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Foreword 


The  objective  of  this  research  is  to  map  lateral  variations  of  regional  wave  Q  in  eastern 
Eurasia  using  data  from  various  sources,  including  the  36  new  CNDSN  stations  in  China. 
The  seismic  waves  utilized  include  high-frequency  regional  phases  Lg,  Pg,  Pn,  and  Sn,  as 
well  as  long-period  Rayleigh  waves.  The  project  was  initiated  by  Dr.  Jiakiang  Xie,  who 
completed  much  of  the  high-frequency  regional  wave  analysis  in  the  first  two  years  of  the 
project.  Following  Dr.  Xie’s  departure  for  a  position  at  AFRL,  the  project  was  continued 
at  FDEO  by  co-PI’s  Dr.  James  Gaherty  and  Dr.  Arthur  Ferner-Fam.  As  a  result,  this 
final  report  consists  of  three  appendices,  each  representing  a  distinct  phase  of  the  project. 
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[i]  We  measure  Q  of  crustal  Lg  waves  in  eastern  Eurasia 
using  the  reliable  two-station  methods.  More  than 
5,000  spectral  ratios  are  collected  over  594  interstation 
paths,  and  are  used  to  determine  values  of  Lg  Q0  and  q  (Q  at 
1  Hz  and  its  frequency  dependence)  over  these  paths.  These 
values  are  used  to  derive  tomographic  models  of  laterally 
varying  Q0  and  r|  with  resolutions  between  4  and  10°.  The 
Qo  model  contains  values  that  vary  between  100  and  900. 
Qo  are  the  lowest  in  the  Tibetan  plateau,  increase  to 
moderate  levels  towards  the  east  and  north,  and  reach 
maxima  in  the  Siberian  and  eastern  Europe  Cratons.  Qo 
values  correlate  well  with  regional  tectonics.  Most  q  values 
are  between  —0.14  and  0.50,  with  a  mean  of  0.18.  These 
values  are  lower  than  those  obtained  previously  using  either 
an  Lg  coda  Q  method,  or  a  method  of  simultaneous 
inversion  of  source  spectra  and  Q.  Citation:  Xie,  J.,  Z.  Wu, 
R.  Liu,  D.  Schaff,  Y.  Liu,  and  J.  Liang  (2006),  Tomographic 
regionalization  of  crustal  Lg  Q  in  eastern  Eurasia,  Geophys.  Res. 
Lett.,  33,  L03315,  doi:10.1029/2005GL024410. 

1.  Introduction 

[2]  The  Lg  is  typically  the  most  prominent  short  period 
seismic  phase  observed  over  continental  paths  at  distances 
greater  than  200  km.  Lg  can  be  treated  as  a  sum  of  higher 
mode  surface  waves,  or  multiple  supercritically  reflected 
S  waves  in  the  crust.  The  quality  factor  (Q)  of  Lg  closely 
resembles  the  crustal  average  of  shear  wave  Q.  The  Lg  Q 
is  known  to  be  highly  variable,  by  up  to  an  order  of 
magnitude,  across  major  continents  [e.g.,  Xie  and  Mitchell, 
1990;  Xie  et  al.,  2004].  Such  large  variations  are  caused  by 
the  fact  that  Q  is  strongly  affected  by  crustal  properties  such 
as  fluid  contents  and  temperature,  both  of  which  can  vary 
significantly  from  one  region  to  another  [ Mitchell  et  al., 
1997].  Reliable  measurements  of  Lg  Q  are  difficult  because 
the  Lg  amplitudes  are  affected  by  complications  of  source 
spectra  and  3D  velocity  structures  whose  effects  are  difficult 
to  fully  account  for,  resulting  in  errors  in  Q  measurement. 
Assuming  that  the  long-duration  coda  of  Lg  is  composed  of 
singly  scattered  Lg  by  crustal  heterogeneities,  one  can 
measure  Lg  coda  Q  using  individual  seismograms  [Xie 
and  Mitchell,  1990;  Mitchell  et  al.,  1997],  The  Lg  coda  Q 
should  grossly  approximate  the  Lg  Q  in  the  areas  sampled 
by  the  scattered  waves. 
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[3]  Among  methods  of  measuring  Q  using  direct  Lg 
waves  [e.g.,  Xie  et  al.,  1996;  Phillips  and  Hartse,  2002], 
the  standard  two-station  method  has  the  advantage  in  that 
it  allows  the  source  spectra  to  be  canceled  in  Q  measure¬ 
ments.  The  reversed  two  station  method  [Chun  et  al., 
1987]  further  allows  the  cancellation  of  site  responses. 
These  methods,  while  being  the  most  reliable,  require  strict 
recording  geometries  and  hence  dense  station  and  event 
coverages. 

[4]  The  eastern  part  of  the  Eurasian  continent  that 
includes  regions  in  and  around  China  has  a  complex 
tectonic  history.  Currently  the  collision  between  the  Indian 
and  Eurasian  plates  to  the  southwest  results  in  the  uplift¬ 
ing  of  the  Tibetan  plateau.  The  extension  along  the  east 
coast  of  China  may  be  related  to  the  underthrusting  of  the 
Pacific  plate  from  the  east.  These  collision  and  under¬ 
thrusting  processes  serve  as  the  primary  driving  forces  of 
the  current  motions  and  deformations  of  a  collage  of 
blocks  with  diverse  histories  of  evolution.  Timing  of  the 
last  significant  tectonic  activity  that  modified  the  crustal 
blocks  ranges  between  Paleozoic  and  current.  Mitchell  et 
al.  [1997]  found  the  Lg  coda  Q  values  throughout  Eurasia 
tend  to  correlate  with  the  time  that  has  elapsed  since  the 
last  major  tectonic  activity.  In  the  past  decades  there  has 
been  a  steady  accumulation  of  digital  Lg  waveforms  from 
eastern  Eurasia  brought  by  the  installation  of  various 
broad  band  seismic  networks  and  the  high  seismicity.  In 
this  paper  we  report  a  tomographic  regionalization  of  Lg 
Q  in  the  region,  obtained  by  applying  the  two  station 
methods  to  more  than  5,000  recently  collected  Lg  spectral 
ratios. 

2.  Data  and  Method  of  Q  Measurement 

[5]  One  hundred  and  sixty  two  broad-band  seismic 
stations  are  used  in  this  study.  The  network  affiliations 
include  the  Global  Seismic  Network  (GSN)  and  various 
national  or  portable  networks  such  as  the  Chinese  Na¬ 
tional  Digital  Seismic  Network  (CNDSN),  the  Kyrghistan 
and  Kazakhstan  seismic  networks,  and  the  three  passive 
networks  deployed  in  the  Tibetan  Plateau  [Xie  et  al., 
2004],  Vertical  component  Lg  waveforms  from  186 
events  between  1988  and  2004  are  retrieved  from  Data 
Management  Centers  of  Incorporated  Research  Institution 
of  Seismology  (IRIS)  and  CNDSN.  Most  events  are 
moderately  sized  (with  magnitudes  clustered  around 
5.5).  Fourier  spectra  of  Lg  are  obtained  using  a  well- 
established  procedure  to  deal  with  noise  and  wave  dis¬ 
persion  [Xie  et  al.,  1996,  2004].  More  than  6,000  Lg 
spectra  are  collected.  From  these,  5,787  pairs  of  spectra 
are  selected  from  two  stations  that  are  (1)  aligned 
approximately  with  the  source,  (2)  separated  far  enough 
(>250  km),  permitting  the  use  of  the  standard  two-station 
method  for  Lg  Q  measurement  [cf.  Xie  et  al.,  2004, 
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Figure  4],  About  a  quarter  of  these  spectral  pairs  further 
satisfy  the  criterion  of  the  reversed  two-station,  two-event 
condition  (i.e.,  the  sources  that  generated  the  ratios  are 
located  on  the  opposite  azimuths  viewed  from  the  stations 
[c.f.  Chun  et  al.,  1987,  Figure  4]),  permitting  estimation 
of  both  inter-station  Q  free  of  site  responses,  and  ratios  of 
site  responses  [Chun  et  al.,  1987].  By  examining  the 
latter  ratios  a  few  stations  are  found  to  have  non-unity 
site  responses  and/or  erroneously  documented  instrument 
responses  at  levels  higher  than  about  2.  These  stations  are 
screened  out  in  the  subsequent  analysis.  This  results  in 
5,265  useful  spectral  ratios  that  collectively  sample  594 
inter-station  paths.  Assuming  that  the  Lg  Q  follow  the 
power-law  frequency  dependence  (Q  =  Qof]  where  Qq  is 
Q  at  1  Hz  and  q  is  the  frequency  dependence),  the  inter¬ 
station  Qq  and  p  can  be  estimated  using  the  following 
concise  relationship  (c.f.  equations  (2)  through  (4)  of  Xie 
et  al.  [2004]): 


In]  — In 

^Av 
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1/2 


.  Mf) 


1/2 


Mf) 


(1 -p)  In/ -In  go,  (1) 


where  VLg  is  the  typical  Lg  group  velocity  of  3.5  km/s,  A,- 
and  Aj  are  the  epicentral  distances  to  the  two  aligned 
stations  where  Lg  spectra  At{f)  and  Aj(f)  are  collected,  A tj- 
is  the  interstation  distance.  In  the  above  equations  the  left 
hand  quantity  is  calculable  from  the  observed  spectral  ratios 
(Ai(f)IAj(f))  and  can  be  averaged  over  repeating  two- 
station  paths  to  provide  stable  estimate  of  Q0  and  p  by  a 
linear  regression.  We  shall  denote  this  quantity  as  “reduced 
spectral  ratio.”  Average  reduced  spectral  ratios  are  calcu¬ 
lated  over  the  594  paths  to  estimate  interstation  Lg  Q0  and 
p.  Typically,  the  lowest  frequencies  that  yield  usable 
reduced  spectral  ratios  are  between  0.1  and  0.3  Hz.  The 
highest  usable  frequencies  are  primarily  controlled  by  the 
signal/noise  ratio  threshold  of  2  used  in  this  study,  and 
typically  vary  between  about  1  and  2.5  Hz.  Figures  1  and  SI 


Figure  1.  Examples  of  average  inter-station  reduced 
spectral  ratios  (circles)  and  the  best  fit  Lg  Q  models  (lines, 
see  equation  (1)).  Station  codes,  interstation  distances  and 
the  estimated  Qq,  p  values  are  written  in  the  panels.  Left 
panel  shows  the  average  ratio  over  an  inter-station  path  in 
southwest  China,  with  sources  located  in  the  same  direction, 
so  the  standard  two-station  method  [e.g.,  Xie  et  al.,  2004]  is 
used.  Right  panel  shows  the  average  ratio  that  samples  an 
inter-station  path  mostly  located  in  eastern  Siberia,  from 
sources  located  in  opposite  directions,  so  the  method  of 
Chun  et  al.  [1987]  is  used. 
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Figure  2.  Map  showing  two-station  paths  used  in  this 
study,  with  colors  indicating  the  values  of  Q0  estimated  for 
the  paths.  An  enlarged  version  of  this  map  is  shown  in 
Figure  S2  in  the  supplemental  material. 


in  the  supplemental  material1  show  examples  of  reduced 
spectral  ratios  and  the  best  fit  Q  models. 

3.  Result  of  Q  Tomography 

[6]  Figures  2  and  S2  show  the  594  two-station  paths  with 
colors  indicating  the  values  of  the  Q0  measurements.  These 
values  tend  to  be  coherent  on  a  large  scale  and  are  input  to  a 
tomographic  back-projection  algorithm  to  invert  for  the 
lateral  variations  in  Qq.  The  back-projection  algorithm  is  a 
well  established,  iterative  algorithm  that  reduces  the  model 
mistfit  (residual)  by  updating  Q0  models  along  one  path  at  a 
time  (for  more  details,  see  Xie  and  Mitchell  [1990]  and  Xie 
et  al.  [2004]).  To  parameterize  the  spatially  varying  Qq,  the 
study  area  is  divided  into  about  2,300  cells  with  constant  Q0 
values  and  a  size  of  2  by  2°.  Figure  3  shows  the  resulting  Q0 
model.  The  errors  and  resolution  associated  with  the  model 
are  estimated  using  the  algorithm  of  Xie  and  Mitchell  [1990] 
and  Xie  et  al.  [2004],  Estimated  random  errors  for  the  Qq 
model  in  Figure  3  are  typically  about  10-15%.  The 
resolution,  as  measured  by  the  point  spreading  functions, 
varies  between  about  4°  in  eastern  China  and  about  10°  in 
higher  latitudes  (>55°N). 

[7]  The  most  striking  low  Qq  region  in  Figure  3  is  that  in 
and  around  the  Tibetan  Plateau  where  Qq  is  between  about 
100  and  200.  Toward  the  east,  Qq  increases  to  about  300  in 
the  Songpan-Ganzi  Belt  and  Qaidam  Basin,  and  to  300-550 
in  the  Yangzi  and  south  China  blocks.  To  the  north  of  these 
regions  there  is  a  band  of  moderate  Qq  regions  (300-450) 
that  covers  the  Tarim  Block,  the  Ordos  and  the  Sino-Korean 
Platforms.  To  the  north  of  this  band  Qq  increases  with 
increasing  latitude  in  the  Altaids,  to  between  about  400 
and  600.  Much  of  the  Kazak  Massif  contains  Qq  values  that 
are  greater  than  600.  Variations  of  Qq  in  the  northernmost 
subregions  can  be  resolved  only  at  scales  of  about  10°  due 
to  sparse  ray  coverage,  with  two  broad  regions  of  high  Qq 
values  between  about  700  and  900  in  the  Siberian  and 
Eastern  Europe  Cratons.  Between  these  regions  lies  the 
Siberian  Trap,  a  province  that  was  affected  by  wide-spread 
volcanism  and  rifting  in  the  late  Paleozoic-Mesozoic  time. 
Values  of  Qq  in  the  Trap  are  between  about  400  and  500, 


'Auxiliary  material  is  available  at  ftp://ftp.agu.org/apend/gl/ 
2005GL024410. 
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Figure  3.  Lg  Q0  map  for  eastern  Eurasia  obtained  in  this  study,  and  simplified  tectonic  boundaries  [e.g.,  Mitchell  et  al., 
1997;  Hearn  et  al.,  2004;  Liang  et  al.,  2004].  Abbreviations  are:  Tibetan  Plateau  (TP),  Tarim  Block  (TB),  Himalayas  (HI), 
Songpan-Ganzi  Belt  (SG),  Southeast  Asia  subplate  (SEA),  Yangtze  block  (YB),  South  China  Block  (SCB),  Sino-Korean 
block  (SK),  Suolun-Xiamulun  block  (SX),  Kazak  Massif  (KM),  Siberia  Craton  (SC),  Siberia  Trap  (ST)  and  Eastern  Europe 
Craton  (EC).  The  Ordos  block  is  centered  at  about  1 17°E  and  37.5°N  [e.g.,  Liang  et  al.,  2004]  and  is  merged  into  the  SK  for 
simplicity.  Note  the  maximum  latitude  plotted  is  5°  less  than  in  Figure  2. 


relatively  lower  than  in  the  Cratons.  In  general,  Lg  Q0 
values  in  Figure  3  exhibit  a  good  correlation  with  the  time 
since  the  last  major  tectonic  events  that  modified  the  blocks, 
as  found  by  Mitchell  et  al.  [1997]  for  Lg  coda  Q0  values. 
The  variations  of  Lg  Q0  in  Figure  3  in  the  southern  and 
eastern  portions  are  more  drastic  than  those  of  Lg  coda  Qo 
shown  in  Plate  1  of  Mitchell  et  al.  [1997].  This  indicates  a 
higher  resolution  of  the  Lg  Q  model  than  the  Lg  coda  Q 
model,  which  is  expected  since  the  Lg  coda  is  composed  of 
scattered  Lg  waves  that  sample  elliptical  areas  rather  than 
great-circles.  Lateral  variations  in  Lg  coda  Q0  is  inherently 
smeared  out  by  this  elliptical  sampling  pattern. 

[8]  The  frequency  dependence  (r|)  values  are  measured 
using  the  slope  of  the  linear  regression  to  the  reduced 
spectral  ratios  (equation  (1);  Figure  1  and  Figure  SI).  In 
theory  fitting  the  slope  (1  —  nq)  is  more  unstable  than  fitting 
the  intercept  (In  Q0)  when  the  data  (reduced  spectral  ratios) 
are  repeatedly  sampled  in  a  narrow  range  of  In  f  near  zero. 
In  practice,  Xie  et  al.  [2004]  state  that  to  obtain  stable 
measurements  of  r|,  a  wide / range  between  about  0.1  and  a 
few  Hz  is  desirable.  Unfortunately,  many  spectral  ratios 
used  in  this  study  are  only  available  in  a  relatively  narrow 
frequency  range  between  about  0.1  and  1  Hz  owing  to  the 
rapid  loss  of  high  frequency  signals  at  the  more  distant 
stations.  We  are  therefore  less  confident  at  the  measured  p 
values  for  many  individual  paths.  We  examine  the  gross 


distribution  of  the  594  r|  measurements  by  calculating  their 
mean  and  standard  deviation,  which  are  0.18  and  0.32, 
respectively.  Therefore  most  p  values  are  between  —0.14 
and  0.50.  The  mean  p  of  0.18  is  lower  than  the  mean  r| 
values  of  0.4-0. 5  obtained  using  Lg  coda  [. Mitchell  et  al., 
1997]  or  a  method  of  simultaneous  inversion  of  source 
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Figure  4.  Map  of  power-law  frequency  dependence  (r|) 
obtained  in  this  study.  The  cell  size  used  in  the  inversion  is 
4  by  4°  and  twice  as  large  as  that  used  in  Qo  inversion 
(Figure  3).  The  larger  size  is  used  to  avoid  spurious 
oscillations  of  r|  caused  by  its  measurement  instability.  An 
enlarged  version  of  this  map  is  shown  in  Figure  S3  in  the 
supplemental  material. 
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spectra  and  Lg  Q  [ Xie  et  al.,  1996].  This  discrepancy  in 
measured  q  values  may  result  from  systematic  bias  in  the 
different  measurement  methods,  each  is  based  on  specific 
assumptions.  The  Lg  coda  Q  method  assumes  a  single, 
frequency  independent  scattering  process.  The  simultaneous 
inversion  of  source  spectra  and  Q  assumes  an  idealized, 
omega-squared  source  model.  The  two-station  method 
assumes  nothing  about  Lg  scattering  or  source  model,  but 
has  a  strict  requirement  of  recording  geometry  which  often 
leads  to  a  large  epicentral  distance  (>1500  km)  of  the  more 
distant  station.  The  assumptions  of  Lg  scattering  or  omega- 
squared  source  model  may  not  be  perfectly  valid,  or  i  |  could 
tend  to  decrease  with  recording  distance  in  addition  to 
varying  laterally  [Xie  et  al,  1996].  Assuming  that  the 
distance  dependence  is  not  the  case  and  r|  values  measured 
in  this  study  are  not  grossly  biased,  we  may  map  the  lateral 
variation  of  q  spatially  which  results  in  broad  regions  with 
low  q  values  (between  —0.15  and  0.25)  throughout  the 
interior  of  the  study  area  (Figures  4  and  S3).  Contrary  to  the 
similarity  in  Q0  models  obtained  using  Lg  (this  study)  and 
Lg  coda  [Mitchell  et  al.,  1997],  the  q  model  in  Figures  4  and 
S3  tends  to  be  significantly  lower  than  those  in  Mitchell  et 
al.  [1997], 

4.  Discussion 

[9]  Lateral  variations  of  Lg  Q  at  1  Flz  obtained  in  this 
study  generally  correlate  with  time  since  the  last  major 
tectonic  event  that  modified  the  crustal  Q  structure.  This 
correlation  is  similar  to  that  found  in  a  previous  study  of  Lg 
coda  Q  [Mitchell  et  al.,  1997].  Values  of  power-law 
frequency  dependence  (p)  of  Lg  Q  measured  in  this  study 
are  systematically  lower  then  those  obtained  previously 
using  Lg  coda,  or  using  a  method  of  simultaneous  inversion 
of  Lg  source  spectra  and  path  Q.  This  discrepancy  may  be 
caused  by  break-downs  of  the  stochastic  Lg  scattering  or 
source  models  assumed  by  the  previous  methods,  or  an 
inherent  distance-dependence  of  p  caused  by  a  distance- 
evolution  of  the  Lg  modal  composition  [Xie  et  al.,  1996], 
Resolving  this  discrepancy  is  important  because  path- 


corrections  at  high  frequencies  are  critically  dependent  on 
accurate  p  values,  and  because  inference  of  layered  crustal 
Q  structures  and  their  causes  can  be  made  using  p  values 
[Mitchell  et  al.,  1997;  Xie  et  al.,  2004],  Future  research 
using  data  with  high  signal/noise  ratios  in  a  wider  frequency 
range  should  enable  us  to  determine  p  values  reliably,  and  to 
explore  the  cause  of  the  p  discrepancy. 
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Abstract  Pn  spectra  are  collected  from  three  PASSCAL  seismic  experiments  in  the 
Tibetan  Plateau  (TP)  over  four  path  groups  to  study  Pn  attenuation.  The  first  and  second 
path  groups  run  southward  from  the  Lop  Nor  Test  Site  (LTS)  to  stations  in  north  and 
south  central  Tibet.  The  third  and  fourth  path  groups  form  two  NNE-oriented  profiles  in 
the  eastern  TP.  Events  recorded  by  the  first  2  path  groups  are  also  recorded  over  central 
Asian  paths  running  westerly  from  the  LTS  to  the  Khyrghistan  network  (KNET).  A  com¬ 
parison  of  Pn  spectra  from  these  events  reveal  that  mantle  lid  under  north  central  TP 
attenuates  Pn  more  severely  than  under  central  Asia.  Apparent  Q{)  and  77  (Pn  Q  at  1  Hz 
and  its  frequency  dependence)  are  estimated  with  a  simplified  geometrical  spreading 
(G.S.T.)  of  A-1'3.  Over  path  group  1  that  heavily  samples  northern  TP,  Q0  and  77  are  esti¬ 
mated  to  be  183±33  and  0.3±0.1.  Over  path  groups  2  and  3  that  sample  either  a  mixture 
of  northern  and  southern  Tibet  or  eastern  Tibet,  Q()  and  77  are  between  about  250-270  and 
0.0-0. 1,  respectively.  Over  the  fourth  path  group  that  sample  the  easternmost  TP  the 
respective  estimates  are  374±51  and  0.3±0.1,  similar  to  the  estimates  of  364  and  0.5  for 
central  Asian  paths  from  LTS  to  KNET.  A  comparison  of  total  Pn  attenuations  that 
include  effects  of  both  G.S.T.  and  Q  in  continental  regions  show  that  they  are  similar  in 
relatively  stable  regions  of  central  Asia,  Scandinavia  and  Canadian  shield.  Within  the  TP, 
Pn  attenuation  is  the  highest  under  north  central  Tibet,  and  decreases  toward  south  and 
east.  In  easternmost  Tibet  Pn  attenuation  becomes  similar  to  those  in  stable  regions.  Lat¬ 
eral  variation  of  Pn  attenuation  inversely  correlates  with  that  of  the  Pn  velocity.  Possible 
causes  of  this  variation  include  (1)  a  thermally  driven  mantle  Qp,  and  (2)  region-specific 
velocity  structures,  which  may  be  characterized  by  different  lid  gradients,  and  density  and 
aspect  ratios  of  3D  scatterers  in  the  vicinity  of  the  Moho. 

Abbreviated  title:  Pn  attenuation  in  Tibet 


Introduction 

Pn  is  the  first  arriving  high-frequency  seismic  phase  at  regional  distances  between 
about  2-15°.  It  is  an  important  phase  to  be  used  in  location,  identification  and  size-estima¬ 
tion  of  seismic  events,  and  in  inferring  the  velocity  and  attenuation  structures  of  the 
Earth’s  uppermost  mantle.  At  close  distances  Pn  is  much  like  a  mantle  head  wave,  which 
evolves  into  mantle  turning  waves  with  increasing  distance  ( e.g Sereno  and  Given, 
1990).  This  peculiar  wave  evolution  causes  the  features  of  Pn,  such  as  its  ray  geometry, 
travel  time  and  amplitude  to  vary  with  distance  in  manners  that  are  highly  dependent  on 
the  velocity  gradient  in  the  mantle  lid.  Simplifications  to  these  variations  are  typically 
introduced  when  Pn  data  is  used  to  study  lid  structures.  Pn  arrival  times  have  been  exten¬ 
sively  used  to  study  the  lateral  variations  of  compressional  wave  velocity  (v p)  in  the  lid 
(e.g.,  Heam  et  al.,  1986).  Most  studies  have  ignored  effects  of  the  mantle  velocity  gradi¬ 
ent  and  treated  Pn  as  a  pure  head  wave  at  all  distances.  A  few  studies,  such  as  Zhao  and 
Xie  (1993)  and  Hearn  et  al.  (2004),  accommodated  the  first-order  effect  of  velocity 


gradient  and  obtained  vp  models  that  vary  both  laterally  and  vertically  in  the  mantle  lid. 

Few  studies  have  been  conducted  to  use  Pn  amplitude  to  estimate  its  attenuation,  or 
Q.  A  fundamental  difficulty  in  such  studies  is  caused  by  the  lack  of  knowledge  of  the 
geometrical  spreading  term  (G.S.T.),  which  describes  the  Pn  amplitude  fall-off  with  dis¬ 
tance  owing  to  the  wavefront  expansion  in  the  mantle  lid.  The  Pn  G.S.T.  is  very  sensitive 
to  details  of  the  lid  velocity  structure.  Sereno  and  Given  (1990),  Xie  (1996),  Lay  et  al. 
(2006)  and  Yang  et  al.  (2007)  estimated  the  G.S.T.  of  Pn  for  Scandinavia,  central  Asia 
and  China  using  synthetic  seismograms  calculated  for  various  ID  velocity  models.  They 
found  even  for  these  simplified  models  the  G.S.T.  varied  with  distance  and  frequency  in 
rather  complex  manners.  In  the  real  3D  Earth  structure,  even  more  complications  are 
expected.  For  example,  scattering  and  multiple  bouncing  (whispering-gallery)  processes 
may  affect  the  Pn  wavetrain  (e.g.,  Menke  and  Richards,  1980,  1983;  Kvaerna  and  Doorn- 
bos,  1991;  Ryberg  et  al.,  1995;  Morozov  and  Smithson,  2000;  Nowack  and  Stacey,  2002; 
Nielsen  and  Thybo,  2003).  To  tackle  the  uncertainty  in  G.S.T.,  different  parameterizations 
have  been  used  in  studies  of  Pn  amplitude  attenuation.  On  the  basis  of  coarse  fits  of  syn¬ 
thetic  seismograms,  Sereno  et  al.  (1988)  and  Xie  and  Patton  (1999)  assumed  a  simplified, 
frequency  independent  G.S.T.  of  A  and  estimated  models  of  Pn  Q  ( Q(f ))  in  Scandi¬ 
navia  and  central  Asia.  They  assumed  a  power-law  frequency  dependence  of  Pn  Q, 

Qif)  =  Qor  ,  (i) 

where  Q0  and  7  are  Pn  Q  at  1  Hz  and  its  frequency  dependence,  respectively.  They 
obtained  similar  Q0  values  of  325  and  364,  and  r]  values  of  ~0.5  for  the  two  regions.  Zhu 
et  al.  (1991)  simultaneously  estimated  a  frequency-dependent  G.S.T.  and  Q(f )  for  Cana¬ 
dian  shield.  Other  authors  such  as  Tinker  and  Wallace  (1997),  who  studied  Pn  attenuation 
at  near-regional  distances  within  about  350  km,  simply  parameterized  Pn  attenuation  by  a 
single  G.S.T.  without  Q. 

The  mantle  lid  under  the  Tibetan  plateau  (TP)  is  known  to  have  highly  laterally 
variable  P  wave  velocities  (vp).  Zhao  and  Xie  (1993),  McNamara  et  al.  (1997)  and  Hearn 
et  al.  (2002)  used  Pn  arrival  time  data  to  invert  for  the  laterally  varying  vp  models.  Their 
v p  models  contain  grossly  similar  features  such  as  a  zone  of  low  vp  under  northern  TP, 
and  high  vp  under  much  of  southern  TP.  Ni  and  Barazangi  (1983)  and  McNamara  et  al. 
(1995)  also  found  that  the  Sn  wave,  which  is  the  shear  wave  that  traverses  the  mantle  lid 
similarly  to  Pn,  seemed  to  suffer  abnormally  high  attenuation  under  northern  TP.  These 
findings  are  used  to  infer  that  the  mantle  lid  is  hot  and  partially  molten  under  a  broad 
region  in  northern  TP,  and  cold  under  southern  TP.  Results  of  Pn  and  Sn  studies  contrib¬ 
uted  to  the  development  of  a  geodynamic  model  in  which  the  Indian  lithosphere  under¬ 
thrusts  in  southern  TP  (Nelson  et  al.,  1996),  and  subducts  into  the  mantle  further  south, 
along  the  latitude  of  ~32°N  (Tilmann  et  al.,  2003). 

Measurements  of  Sn  attenuation  under  TP  is  limited  to  a  qualitative  comparison  of 
the  relative  amplitude  of  Sn  to  that  of  P  coda.  By  contrast,  the  amplitude  spectra  of  the 
Pn  wavetrain  can  be  quantitatively  calculated  at  the  distances  where  Pn  is  the  first  arrival. 
In  this  paper  we  present  analyses  of  the  attenuation  of  Pn  spectra  recorded  during  three 
PASSCAL  experiments  in  TP.  We  first  analyze  Pn  spectra  from  three  events  that 
occurred  in  1994  or  1999  in  the  Lop  Nop  Nor  Test  Site  (LTS)  of  eastern  Tarim  Basin, 
which  are  simultaneously  recorded  by  the  INDEPTH  (International  Deep  Profiling  of 
Tibet  and  the  Himalaya)  II  or  III  stations  to  the  south  crossing  the  TP,  and  Khyrghistan 
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network  (KNET)  to  the  west  crossing  eastern  Tienshan.  A  comparison  of  raw  Pn  spectra 
recorded  by  the  INDEPTH  and  KNET  stations  provides  direct  evidence  that  the  mantle 
lid  under  the  TP  is  more  attenuative  than  that  under  eastern  Tienshan.  We  then  quantify 
the  Pn  attenuation  along  paths  to  INDEPTH  networks  by  estimating  apparent  Pn  Q  under 
a  simplified  G.S.T.  and  omega-squared  source  models.  This  is  followed  by  two-station 
measurements  of  apparent  Pn  Q  along  two  PASSCAL  profiles  in  eastern  Tibet,  deployed 
during  the  1991-1992  Tibetan  Plateau  experiment.  We  report  a  spatial  trend  for  the 
apparent  Pn  Q  at  1  Hz  ( Q0 )  to  be  low  in  north  central  TP,  and  to  increase  both  southward 
and  eastward.  This  trend  correlates  with  the  spatial  variation  of  Pn  velocity  measured  in 
previous  studies.  Finally  we  compare  the  total  Pn  attenuations  reported  for  various  conti¬ 
nental  regions  in  the  world,  by  taking  into  account  both  the  G.S.T.  and  Q  model.  We  find 
that  in  the  Canadian  and  Scandinavian  shields,  central  Asia,  and  easternmost  TP  where 
Pn  velocity  is  high  (>  8.0  km/s),  the  total  Pn  attenuations  are  similar  to  one  another  and 
much  weaker  than  that  under  the  northern  TP.  We  infer  possible  causes  of  the  regional 
variation  of  Pn  attenuation,  including  variations  in  either  a  temperature- sensitive  lid  Qp 
and/or  in  fine-scale  velocity  structures  in  the  vicinity  of  the  Moho. 

Data  Processing 

Figure  1  shows  all  seismic  stations  in  Tibet  deployed  during  the  three  PASSCAL  experi¬ 
ments.  The  1991-1992  Sino-US  Tibetan  Plateau  experiment  consisted  of  11  broad-band 
stations  that  were  deployed  between  July,  1991  and  July,  1992  (Owens  et  al.,  1993)  in  the 
eastern  part  of  the  plateau.  The  INDEPTH  II  experiment  was  conducted  in  southern 
Tibet  between  May  and  October,  1994  (Nelson  et  al.,  1996)  with  15  broad-band  or  short- 
period  stations  deployed  along  an  NNE-oriented  profile.  The  INDEPTH  III  experiment 
was  conducted  in  central  Tibet  between  July,  1998  and  June,  1999.  Forty-seven  broad  or 
intermediate  band  and  15  short-period  stations  were  deployed  along  a  profile  across  the 
Banggong-Nujiang  Suture  (BNS)  in  a  NNW  direction  (Figure  1;  also  c.f.  Rapine  et  al., 
2003).  Stations  in  all  three  experiments  are  three-component,  with  sampling  rates  of  20 
s-1  or  higher.  Figure  1  also  shows  regional  events  used  in  this  study.  Figure  2  shows  an 
example  of  the  record  sections  containing  Pn. 

To  obtain  stable  estimates  of  Pn  amplitude  spectra,  a  finite-length  wavetrain  that 
included  Pn  and  its  coda  was  previously  used  by  Sereno  et  al.  (1988),  Zhu  et  al.  (1991) 
and  Xie  and  Patton  (1999).  In  this  study  we  follow  the  procedure  by  these  authors  to 
compute  average  spectra  of  Pn  and  Pn  coda,  using  spectra  from  a  series  of  overlapping 
windows  applied  to  the  vertical-component  seismograms.  These  windows  have  a  constant 
length  of  4.5  s,  a  20%  taper  and  50%  overlap.  Only  spectra  with  signal  to  noise  ratios 
larger  than  2  are  used.  The  maximum  number  of  the  windows  for  a  given  seismogram  is 
determined  by  the  following  criteria:  (a)  the  total  length  of  Pn  wavetrain  sampled  does 
not  exceed  12.5  s  (i.e.,  the  maximum  window  number  is  5),  (b)  at  closer  distances  (A  < 
1000  km)  the  last  window  does  not  include  signals  with  group  velocities  slower  than  6.6 
km/s  so  that  the  crustal  Pg  wave  is  excluded;  and  (c)  at  large  distances  (A-1500  km)  the 
last  window  does  not  include  the  possible  reflections  from  the  440  km  and  660  km  dis¬ 
continuities,  previously  observed  in  Tibet  (8301  Program  Group,  1988).  Figure  3  shows 
an  example  of  the  windows  used,  and  the  resulting  individual  and  average  Fourier  spectra 
of  the  Pn  wavetrain.  The  effect  of  ambient  noise  is  reduced  by  subtracting  the  signal 
power  by  the  noise  power.  Instrument  responses  are  then  removed.  We  shall  refer  to  the 
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average  spectra  thus  calculated  for  the  Pn  wavetrain  as  "Pn  spectra"  in  the  rest  of  the 
paper. 


Features  of  Pn  Spectra  from  Events  in  the  Lop  Nor  Test  Site 

Spectra  from  the  1994  Explosions 

During  the  INDEPTH  II  experiment,  two  underground  nuclear  explosions  occurred 
in  the  LTS  of  eastern  Tarim  Basin  (Figure  1).  Pn  wavetrains  were  recorded  by  all  12 
operational  INDEPTH  stations  for  the  October  7  (mb= 5.9)  explosion  (Figure  2),  and  by  7 
stations  for  the  June  10  (mb= 5.7)  explosion.  Pn  spectra  from  both  explosions  were  also 
recorded  by  the  Khyrghistan  network  (KNET)  to  the  west  (Figure  1),  and  studied  by  Xie 
and  Patton  (1999).  The  paths  to  the  KNET  stations  traverse  to  the  eastern  Tienshan.  In 
comparison,  the  paths  to  the  INDEPTH  II  network  run  southward  and  sample  the  TP. 
Both  groups  of  paths  have  similar,  far-regional  distances  between  about  1050  and  1550 
km  (Figure  1).  Any  anisotropic  Pn  radiation  patterns  are  expected  to  be  minimal.  There¬ 
fore  we  may  directly  compare  the  Pn  spectra  to  see  the  gross  attenuation  effect  of  the 
mantle  lid  under  Tienshan  and  the  Tibetan  Plateau.  The  top  two  rows  of  Figure  4  show 
the  Pn  ground  displacement  spectra  from  both  explosions,  scaled  by  A  to  reduce  the 
effects  of  G.S.T. 

A  prominent  feature  of  the  Pn  spectral  amplitudes  from  both  explosions  is  their 
variability  within  the  networks,  particularly  within  the  KNET.  Both  networks  have  aper¬ 
tures  of  less  than  about  350  km  compared  to  the  average  path  lengths  of  >  1050  km.  Yet 
the  sample  standard  deviation  is  typically  at  the  level  of  -60%  of  the  averaged  spectra  for 
the  KNET  stations,  and  range  between  -60%  (at  frequencies  <  1  Hz)  and  -30%  (at  fre¬ 
quencies  >  1  Hz)  for  the  INDEPTH  stations.  For  both  networks  the  amplitude  variations 
are  not  monotonic  with  distance,  so  path  attenuation  from  finite  Q  is  not  the  primary 
cause  of  these  variations.  Xie  (1996)  and  Xie  and  Patton  (1999)  concluded  that  the  most 
likely  cause  for  the  drastic  Pn  amplitude  variations  across  the  KNET  is  a  laterally  varying 
Moho  topography  that  likely  occurs  along  the  Tienshan  range.  Xie  (1996,  2003)  con¬ 
ducted  ray  tracings  in  2D  velocity  models  for  Central  Asia  in  and  around  the  Tienshan 
region  to  demonstrate  the  potential  effect  of  Moho  topography.  Significant  amplitude 
variations  of  Pn  over  short  distance  ranges  seem  to  be  common  in  other  mountainous 
regions  in  the  world  (e.g.,  Tinker  and  Wallace,  1997).  The  use  of  a  finite  duration  wave 
train  to  calculate  Pn  spectra  has  not  sufficiently  suppressed  its  variability.  Therefore  to 
study  the  gross  effect  of  whole-path  attenuations,  we  take  the  network  average  of  Pn 
spectra,  as  shown  in  the  right  columns  of  Figure  4. 

The  network- averaged  Pn  spectra  from  the  well  recorded,  October  7  explosion  (top 
right  panel  of  Figure  4)  start  at  similar  levels  at  similar  distances  for  KNET  and 
INDEPTH  II  at  the  lowest  frequency  (-0.4  Hz).  With  increasing  frequency  the  average 
spectrum  from  the  INDEPTH  II  network  rapidly  falls  below  that  from  the  KNET.  At  -1 
Hz  the  separation  between  the  two  spectra  reaches  a  factor  of  4.  At  higher  frequencies  the 
separation  increases  monotonically  with  frequency  and  overwhelms  the  variabilities 
within  the  networks  as  measured  by  the  standard  deviations.  At  4  Hz  (the  highest  fre¬ 
quency  available  for  INDEPTH  spectra),  the  separation  reaches  a  factor  of  -35.  Network 
averaged  spectra  from  the  June  10  explosion  (right  panel  of  row  2  of  Figure  4)  have  a 
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similar,  frequency-dependent  separation,  although  the  smaller  event  size  causes  both 
fewer  recordings  and  narrower  recording  bands.  Since  the  recording  distances  to  KNET 
and  INDEPTH  networks  are  similar,  the  frequency  dependent  separations  between  the 
network- averaged  Pn  spectra  from  both  explosions  mean  that  Pn  attenuation  under  the 
Tibetan  plateau  and  Tienshan  is  similar  at  low  frequencies  (~0.4  Hz),  but  becomes 
increasingly  different  with  increasing  frequency. 

Spectra  from  the  1999  Earthquake 


Event  "99.030"  in  Figure  1  is  an  mb= 5.9  earthquake  on  January  30,  1999,  which  is 
located  nearby  the  1994  explosions.  The  event  is  simultaneously  recorded  by  10  KNET 
and  38  INDEPTH  III  stations.  Of  the  latter,  35  are  deployed  along  a  main  profile  that 
crosses  the  BNS,  and  3  are  deployed  on  the  BNS  but  outside  the  main  profile  (Figures  1 
and  5).  Pn  velocity  analyses  reveal  that  the  BNS  is  an  approximate  boundary  separating 
regions  of  abnormally  slow  and  fast  mantle  lid  velocities  (Zhao  and  Xie,  1993;  McNa¬ 
mara  et  al.,  1997;  Hearn  et  al.,  2002;  Xie,  2003).  So  in  this  study  we  have  divided  the 
INDEPTH  III  stations  into  two  groups.  Group  1  (denoted  as  III- 1 )  contains  15  stations 
south  of  the  BNS  and  group  2  (denoted  as  III-2)  contains  23  stations  north  of,  and  on,  the 
BNS  (Figure  5).  Rows  3  and  4  of  Figure  4  show  the  Pn  spectra  from  KNET  and  the  two 
groups  of  INDEPTH  III  stations.  As  in  the  case  of  Pn  spectra  from  the  explosions  (top 
two  rows  in  Figure  4),  the  spectra  from  the  earthquake  (bottom  two  rows)  also  exhibit 
substantial  variations  within  each  network.  The  averages  of  Pn  spectra  from  KNET  and 
INDEPTH  III  stations  also  exhibit  a  frequency-dependent  separation,  indicating  more 
severe  attenuation  of  the  Tibetan  mantle  lid  at  higher  frequencies.  There  are  also  notable 
differences  between  explosion  and  earthquake  spectra.  Whereas  the  explosion  spectra 
contain  spectral  overshoot  effects  near  the  source  corner  frequencies  (e.g.,  above  1  Hz  for 
the  October  7  explosion),  the  earthquake  spectra  (bottom  rows)  exhibit  no  overshoot  and 
decay  more  gently  (see  Xie  and  Patton,  1999,  for  details). 

Stochastic  Modeling  of  Pn  Spectra 


To  quantitatively  estimate  Pn  Q  and  source  spectra,  we  use  the  stochastic  modeling 
of  Pn  spectra  used  by  Sereno  et  al.  (1988)  and  adapted  by  Xie  and  Patton  (1999).  For  con¬ 
venience  we  briefly  summarize  the  modeling  here.  We  assume  that  Pn  spectra,  A(/),  can 
be  modeled  by 


A(f)  =  S(f)R(0)G(  A)  exp 


K 


f  A 


VQ(f ) 


m) 


(2) 


where  6,  A,  V  are  the  azimuth,  distance  and  Pn  group  velocity,  respectively.  Q(f)  is  Pn  Q 
(equation  (1)).  R{6 )  is  the  source  radiation  pattern.  X(f)  is  an  error  term  that  represents 
systematic  errors  such  as  a  non-unity  site  response,  and  random  errors  caused  by  ampli¬ 
tude  fluctuation.  S(f )  is  the  Pn  source  spectrum,  which  is  given  by  the  Brune’s  model 
and  Modified  Mueller-Murphy  (MMM)  model  for  earthquakes  and  explosions,  respec¬ 
tively  (equation  (2)  of  Xie  and  Patton,  1999).  Both  models  have  f~ 2  asymptotic  decays 
at  high-frequencies,  and  are  grossly  characterized  by  a  seismic  moment  (M0)  and  corner 
frequency  (fc).  The  MMM  model  for  explosions  also  has  a  spectral  overshoot  controlled 
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by  a  parameter  /?,  which  is  set  as  1.0  by  Xie  and  Patton  (1999).  In  equation  (2),  G(A)  is 
the  geometrical  spreading  term  (G.S.T.)  and  takes  the  form, 

G(A)  =  Ao1(A0/A)M  (3) 

where  A0  is  a  reference  distance,  and  m  is  the  decay  rate  of  A(f)  at  large  distances 
(A  >  A0).  As  mentioned  in  the  Introduction,  a  well  known  problem  in  modeling  Pn  is  that 
the  values  of  A()  and  m  are  uncertain  because  of  the  unknown  details  of  the  lid  velocity 
structure.  Based  on  theoretical  and  empirical  considerations,  Sereno  et  al.,  (1988)  used 
A0  and  m  of  1  km  and  1.3,  respectively,  for  Scandinavia.  Xie  (1996)  found  that  these  val¬ 
ues  were  also  within  the  ranges  obtained  using  synthetic  Pn  waveforms  in  various  ID 
velocity  models  for  central  Asia.  Therefore  Xie  and  Patton  (1999)  used  these  A0  and  m 
values  to  study  Pn  spectra  from  several  nearly  colocated  explosions  in  the  LTS,  recorded 
by  the  KNET  in  the  west,  and  the  Kazakhstan  network  (KZNET)  in  the  north.  To  obtain 
robust  estimates  of  source  and  path  parameters  with  a  minimum  trade-off  among  them, 
Xie  and  Patton  first  estimated  source  fc  values  using  the  empirical  Green’s  function 
method,  and  path-variable  Pn  Q  using  a  two-station  and  a  composite  event  method.  These 
estimates  were  then  input  as  a  priori  knowledge  to  an  Bayesian  inverse  method  to  mea¬ 
sure  M0  and  refined  fc  values  for  all  explosions.  The  Pn  Q  values  by  Xie  and  Patton 
(1999)  are  apparent  because  of  the  specific  G.S.T  used,  but  are  robust  because  the  multi¬ 
ple-event,  multiple-constraints  used.  These  values  will  be  used  as  a  priori  knowledge  in 
this  study. 


Estimates  of  Apparent  Pn  Q  to  Indepth  Stations 

The  October  7  explosion  was  recorded  by  8  KNET  stations.  The  average  Q0  and  ij 
to  these  stations  were  estimated  to  be  381  and  0.4,  respectively  (Xie  and  Patton,  1999). 
We  input  these  values  to  the  Bayesian  method  to  estimate  source  M0,  fc  and  path-vari- 
able  <20,  7  values  to  the  INDEPTH  II  stations  (Equation  (2),  also  see  equation  (13)  of  Xie 
and  Patton,  1999).  Figure  6  (row  1)  shows  the  fit  of  the  resulting  source  spectral  and  path 
Q  models  to  the  observed  spectra  from  the  INDEPTH  II  network  and  KNET.  The  esti¬ 
mated  M0,  fc  values  are  the  same  as  those  previously  obtained  by  Xie  and  Patton  (1999; 
3.2  x  1015  Nm  and  2.6  Hz,  respectively).  The  estimated  average  Q0  value  to  the 
INDEPTH  II  stations  is  253±53,  which  is  30%  lower  than  the  value  of  364  to  all  12 
KNET  stations  estimated  by  Xie  and  Patton  (1999).  The  average  q  to  INDEPTH  II  sta¬ 
tions  is  0.0±0.1,  which  is  much  lower  than  the  average  q  of  0.5  estimated  to  the  KNET 
stations.  The  June  10  explosion  is  relatively  poorly  recorded  as  mentioned  in  a  previous 
section.  We  are  not  confident  at  estimating  Q  using  that  event. 

Earthquake  99.030  was  recorded  by  all  12  KNET  stations  and  38  INDEPTH  III  sta¬ 
tions.  We  assume  that  for  this  earthquake  the  source  radiation  pattern  (R(6 >)  in  equation 
(2))  does  not  vary  significantly  to  the  former  and  latter  stations.  This  assumption  is  justi¬ 
fied  quantitatively  in  the  Appendix.  The  assumption  allows  us  to  use  the  average  Q0,  q 
values  to  KNET  estimated  by  Xie  and  Patton  (1999;  364  and  0.5,  respectively)  as  a  priori 
information  in  the  Bayesian  inversion  in  the  same  way  as  we  did  for  explosions.  The 
resulting  estimates  of  source  M{]  and  fc  values  are  2.49  x  1016  Nm  and  0.8  Hz,  respec¬ 
tively.  Our  fc  estimate  has  an  uncertainty  of  about  0.1  Hz  and  is  very  consistent  with  the 
estimated  fc  range  between  0.72  and  0.78  Hz  for  this  event,  obtained  separately  using  an 
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empirical  Green’s  function  method  by  Dr.  Mark  Fisk  (M.  Fisk,  written  communication, 
September,  2006).  The  average  Q0  and  77  values  to  the  southern  station  group,  III- 1 ,  are 
estimated  to  be  255±48  and  0.1±0.1.  These  values  are  similar  to  those  to  the  INDEPTH  II 
stations  estimated  above.  By  contrast,  the  average  Q0  and  77  values  to  station  group  III-2, 
located  to  the  north  of  BNS,  are  estimated  to  be  183±33  and  0.3±0.1,  respectively. 

The  Qo  and  77  estimates  presented  in  this  section  are  those  along  long  paths  from 
LTS  and  contain  small  crustal  segments  in  both  the  source  and  receiver  end,  and  short 
mantle  segments  in  the  Tarim  Basin.  The  two-station  method  (Xie  et  al.,  2004),  when 
applicable,  can  minimize  effects  of  both  attenuation  outside  the  networks  and  the  source 
radiation.  We  attempted  to  use  that  method  to  estimate  Pn  Q  along  the  INDEPTH  pro¬ 
files  but  obtained  very  unstable  (negative)  Q  values.  This  is  caused  by  the  large  variability 
of  Pn  spectral  amplitudes  within  the  INDEPTH  II  and  III  networks  (Figure  4)  which 
means  errors  (. X(f )  in  equation  (2))  are  large.  Xie  et  al.  (2004)  pointed  out  that  at  a  given 
frequency,  the  effect  of  error  ( X(f ))  in  Q  estimates  could  be  suppressed  by  using  large 
measuring  distances.  Unfortunately  the  INDEPTH  profiles  are  short  (<  350  or  400  km; 
see  Figure  1),  resulting  in  small  inter-station  distances  that  are  not  adequate  to  suppress 
the  large  X(f). 

Apparent  Pn  Q  along  Passcal  Profiles  in  Eastern  Tibet 

Stations  deployed  during  the  1991-1992  Sino-US  Tibetan  Plateau  experiment 
(Owens  et  al.,  1993)  recorded  Pn  spectra  from  four  regional  earthquakes  (Figures  1  and 
5)  that  are  approximately  in  line  with  profiles  SANG-TUNL  and  SANG/GANZ-MAQI  in 
eastern  Tibet.  Inter-station  distances  along  these  profiles  are  generally  large  so  the  effect 
of  the  error  term  (X(  f))  in  two-station  Q  measurements  should  be  reduced  (last  section). 
We  use  Pn  spectra  along  these  profiles  to  estimate  spectral  ratios  from  many  pairs  of  two- 
stations.  These  ratios  are  then  averaged  along  profiles  SANG-TUNL,  SANG/GANZ- 
MAQI  and  WNDO-TUNL  to  estimate  Q{)  and  77  using  equation  (4)  of  Xie  et  al.  (2004), 
replacing  V Lg  by  a  Pn  phase  velocity  of  8.0  km/s.  Profile  SANG/GANZ-MAQI  is  merged 
from  two  subprofiles  SANG-MAQI  and  GANZ-MAQI.  Profile  WNDO-TUNL  is  actu¬ 
ally  the  northern  half  of  profile  SANG-TUNL;  it  is  used  to  see  if  Pn  Q  changes  from 
south  to  north  along  profile  SANG-TUNL  that  crosses  the  BNS.  Figure  7  shows  the 
stacked  spectral  ratios  (SSRs)  along  the  three  profiles  and  the  best  fitting  models  of  Q0 
and  77.  At  most  frequencies  (>  0.7  Hz)  the  SSRs  closely  follow  the  straight  lines  repre¬ 
senting  the  best  fitting  models.  The  SSRs  at  frequencies  lower  than  0.7  Hz  are  abnormally 
small  along  all  profiles.  These  may  be  caused  by  a  trend  of  increasing  effects  of  modeling 
errors  (X(f))  in  Q  estimates  toward  lower  frequencies,  as  predicted  by  Xie  et  al.  (2004). 
The  <2o  and  7  values  estimated  along  profile  SANG-TUNL  are  278±23  and  0.1  ±0.1, 
respectively  (Figure  7),  similar  to  those  estimated  for  paths  from  Lop  Nor  to  the 
INDEPTH  stations  in  southern  Tibet  (last  section).  There  is  no  significant  north-south 
change  of  Q0  and  77  values  along  profile  SANG-TUNL  since  the  values  estimated  along 
its  northern  half  (profile  WNDO-TUNL)  are  273±32  and  0.0±0.1  (Figure  7),  virtually  the 
same  as  those  estimated  along  the  entire  profile.  The  Q0  and  77  estimated  along  profile 
SANG/GANZ-MAQI,  which  is  located  in  the  easternmost  TB,  are  about  374±57  and 
0.3±0.1,  similar  to  the  values  of  364  and  0.5  estimated  for  paths  from  the  Lop  Nor  to 
KNET  (Xie  and  Patton,  1999,  also  see  the  last  section). 
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A  Comparison  of  Pn  Attenuation  in  Continental  Regions 

Figure  5  summarizes  all  apparent  Q0  and  rj  values  estimated  in  this  study  and  Xie 
and  Patton  (1999)  for  regions  in  and  around  the  Tibetan  plateau  and  Tienshan.  Paths  from 
the  LTS  to  north  central  Tibet  are  associated  with  the  lowest  Q0  (-183).  Higher  Q0  val¬ 
ues  are  found  along  paths  from  LTS  to  south  central  Tibet,  and  along  the  two  profiles  in 
eastern  Tibet.  These  Q0  and  q  values  allow  us  to  compare  Pn  attenuation  in  the  TP  to 
those  in  other  continental  regions.  As  mentioned  in  the  Introduction,  there  have  only 
been  few  previous  publications  documenting  Pn  attenuation  on  continents  at  far-regional 
distances.  The  published  Pn  Q  values  are  obtained  using  different  G.S.T.  (Sereno  et  al., 
1988;  Zhu  et  al.,  1991;  Xie  and  Patton,  1999).  To  compare  them  with  those  obtained  in 
this  study,  we  calculate  models  of  total  Pn  attenuation  in  the  frequency  domain  at  a  refer¬ 
ence  distance  of  1200  km  (Figure  8).  The  calculation  simultaneously  accounts  for  the 
G.S.T.  used,  and  the  estimated  Pn  Q.  Pn  attenuation  models  are  similar  for  the  relatively 
stable  regions  where  mantle  Pn  velocities  are  higher  than  8.0  km/s,  including  the  Cana¬ 
dian  shield,  Scandinavia,  central  Asia  and  easternmost  Tibet  (omitted  in  Figure  8;  see 
caption).  This  is  so  despite  that  the  G.S.T.  used  for  the  Canadian  Shield  by  Zhu  et  al. 
(1991)  has  a  frequency-dependent  form  of  m=l.  07  +  0. 035/,  rather  than  m=1.3  used  for 
other  regions. 

The  lowest  model  (corresponding  to  the  strongest  Pn  attenuation)  in  Figure  8  is  that 
for  paths  from  LTS  to  northern  Tibet  (model  "N.  TP").  At  all  frequencies  the  model  is 
below  those  for  the  stable  regions,  but  the  differences  are  more  drastic  toward  higher  fre¬ 
quencies.  At  4  Hz  the  differences  are  more  than  a  factor  of  100.  The  second  lowest 
model  in  Figure  8  is  Model  C.  TP.  which  at  most  frequencies  is  much  below  those  for  sta¬ 
ble  regions,  while  being  somewhat  higher  than  model  N.  TP.  There  are  two  possible 
causes  for  models  N.  T.P  and  C.  T.P  to  be  lower  than  models  for  stable  regions.  The  first 
one  is  that  the  intrinsic  Qp  under  northern  Tibet  may  be  low.  This  cause  is  empirically 
supported  by  the  monotonic  increase  of  the  differences  between  the  Tibetan  and  other 
models  with  frequency  in  Figure  8.  Such  frequency  content  shift  among  seismic  spectra  is 
a  typical  indicator  of  a  varying  Q  (e.g.,  Ruzaikin  et  al.,  1977;  Fan  and  Lay,  2002;  Xie  et 
al.,  2002).  A  physical  reasoning  for  this  cause  is  based  on  the  fact  that  intrinsic  Qp  is 
known  to  be  strongly  dependent  on  temperature  (e.g.,  Mitchell,  1995).  As  mentioned  in 
the  Introduction,  analyses  of  lateral  variations  in  Pn  velocity  have  led  to  the  inference  that 
the  temperature  in  the  mantle  lid  under  north  central  Tibet  is  abnormally  high,  likely 
reaching  the  solidus  and  causing  partial  melting.  By  contrast,  the  lid  temperature  beneath 
southern  Tibet  (south  of  the  BNS  in  Figures  1  and  5)  is  abnormally  low  owing  to  the 
underthrusting  Indian  lithosphere.  Since  model  N.  TP  most  heavily  samples  the  hot  lid 
under  north  central  Tibet,  it  is  expected  to  be  most  affected  by  the  low  Qp  caused  by  high 
temperature  there.  Model  C.  T.P.  is  developed  for  paths  from  LTS  to  southern  Tibet  (Fig¬ 
ure  5)  and  samples  a  mixture  of  the  hot  lid  under  north  central  Tibet  and  cold  lid  under 
southern  Tibet,  so  the  average  Qp  is  expected  to  be  higher  than  that  of  model  N.  T.P, 
making  model  C.  T.P.  somewhat  above  N.  T.P.  These  expected  Qp  levels  well  explain 
why  models  N.T.P.  and  C.T.P.  are  the  first  and  second  lowest  in  Figure  8.  The  second 
possible  cause  for  the  T.P.  models  to  be  low  is  a  regional  variation  of  the  near  Moho- 
velocity  structure.  Studies  using  synthetic  seismograms  have  shown  that  the  Pn 
wavetrain  may  be  composed  of  diving  and  whispering-gallery  waves  whose  propagation 
efficiency,  or  G.S.T.,  may  be  sensitive  to  details  of  the  velocity  structure  near  the  Moho 
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(Menke  and  Richards,  1980;  1983;  Sereno  and  Given,  1990;  Kvaema  and  Doornbos, 
1991).  In  particular,  a  positive  lid  velocity  gradient  may  retain  more  Pn  energy  by  refrac¬ 
tion.  Likewise  an  enhanced  distribution  of  horizontally  elongated  scatterers  in  the  mantle 
lid  or  lower  crust  may  strengthen  Pn  by  forward- scattering  (e.g.,  Ryberg  et  al.,  1995, 
Morozov  and  Smithson,  2000;  Nowack  and  Stacey,  2002;  Neilsen  and  Thybo,  2003).  As 
will  be  elaborated  below,  we  can  not  rule  out  this  second  cause  because  we  fundamentally 
lack  knowledge  of  detailed  regional  velocity  structures  involved. 

Discussion 

In  this  study  we  have  measured  apparent  Pn  Q  with  a  simplified  geometrical 
spreading  term  (G.S.T.)  of  A-1'3.  That  G.S.T.  was  used  in  the  past  in  Scandinavia  and  cen¬ 
tral  Asia  (Sereno  et  al.,  1990;  Xie  and  Patton,  1999)  and  is  within  a  range  of  G.S.T. 
obtained  using  synthetic  Pn  waveforms  in  various  ID  velocity  models  for  central  Asia 
(Xie,  1996).  We  can  not  improve  that  G.S.T.  using  synthetic  seismogram  calculations 
because  we  do  not  have  detailed  knowledge  of  the  velocity  structure  in  the  vicinity  of  the 
Moho.  For  example,  previous  Pn  velocity  models  for  the  T.P  lack  the  resolution  for 
inferring  whether  there  is  an  abnormally  negative  sub-Moho  gradient  under  north-central 
Tibet  (e.g.,  8301  Group,  1988;  Zhao  and  Xie,  1993;  McNamara  et  al.,  1997).  Inferring 
features  of  3D  scatterers,  such  as  their  density  and  aspect  ratios  that  may  strongly  affect 
Pn  G.S.T.  (end  of  last  section)  would  be  even  more  difficult  as  it  requires  even  higher  res¬ 
olution.  Zhu  et  al.  (1991)  proposed  to  use  the  observed  Pn  spectral  ratios  to  simultane¬ 
ously  fit  a  frequency-dependent  G.S.T.  and  Pn  Q  model.  A  concern  about  that  approach 
is  that  there  is  likely  a  trade-off  between  the  estimated  G.S.T.  and  Q  (Sereno  et  al.,  1990; 
Xie  and  Patton,  1999).  To  see  if  such  trade-off  exists  in  our  dataset,  we  conduct  numeri¬ 
cal  tests  in  which  we  fit  two-station  spectral  ratios  along  the  Passcal  profiles  in  eastern 
Tibet  by  various  combinations  of  G.S.T.  and  Q,  and  examine  how  much  the  resulting 
residuals  vary.  Figure  9  shows  an  example  of  such  tests,  in  which  stacked  spectral  ratios 
along  profile  SANG-TUNL  are  fit  by  the  G.S.T.  of  m  =  1.07  +  0.0035/  estimated  by 
Zhu  et  al.  (1991),  together  with  a  Q0  and  r/  of  241  and  0.20.  The  residual  of  this  new  fit  is 
0.6646.  In  comparison,  the  residual  obtained  earlier  using  m  =  1.  3,  <2o  =  278  and  tj  =  0.1 
(top  panel  of  Figure  7;  replotted  in  Figure  9)  is  0.6660.  These  two  residuals  differ  by  only 
0.2%,  suggesting  a  nearly  complete  trade-off  between  G.S.T.  (m)  and  Q  ( Q0  and  rj).  Since 
we  will  unlikely  be  able  to  acquire  precise  knowledge  of  the  Pn  G.S.T.  and  its  regional 
variations  in  the  near  future,  the  estimated  Pn  Q  will  be  apparent  Q  which  only  approxi¬ 
mates  the  true  Q.  Such  apparent  Q  does  not  provide  input  for  a  precise  estimate  of  lid 
temperature  using  a  Q-T  relationship  (e.g.,  Mitchell,  1995).  In  the  practical  path  correc¬ 
tions,  total  Pn  attenuation  that  combines  the  effects  of  Pn  Q  and  the  respective  G.S.T. 
should  be  used. 


Conclusions 

Pn  attenuation  is  studied  using  data  from  three  PASSCAL  experiments  in  the 
Tibetan  Plateau  (TP),  including  the  1991-1992  Sino-US  TP,  the  INDEPTH  (International 
Deep  Profiling  of  Tibet  and  the  Himalaya)  II  and  III  experiments.  Three  seismic  events 
in  the  Lop  Nor  Test  Site  (LTS)  are  simultaneously  recorded  by  the  INDEPTH  II  and  III 
networks  to  the  south  and  the  Khyrghistan  network  (KNET)  to  the  west,  at  similar  far- 
regional  distances  (>  1050  km).  A  comparison  of  the  recorded  Pn  spectra  yields  direct 
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evidence  that  the  mantle  lid  under  the  TB  is  more  attenuative  than  that  under  central  Asia 
at  high  frequencies  (>  1  Hz).  At  4  Hz,  the  difference  in  attenuations  cause  the  INDEPTH 
spectra  to  be  about  a  factor  of  35  smaller  than  the  KNET  spectra. 

Under  a  simplified  Pn  geometrical  spreading  of  A-1'3,  we  measure  Q{)  and  i)  (appar¬ 
ent  Pn  Q  at  1  Hz  and  its  frequency  dependence,  respectively)  along  various  path  groups. 
For  paths  from  LTS  to  northern  Tibet  we  estimate  Q0,  tj  to  be  1 83+33  and  0.3±0.1, 
respectively.  For  paths  to  southern  Tibet,  and  along  a  profile  (SANG-TUNL)  in  eastern 
Tibet  deployed  during  the  1991-1992  experiment,  Q0  increases  to  above  250,  while  77 
decreases  to  about  0.0-0. 1.  Along  a  profile  in  the  easternmost  Tibet  (SANG/GANZ- 
MAQI),  <2o  increases  further  to  374±51,  while  77  is  0.3±0.1.  These  values  are  similar  to 
those  of  364  and  0.5  measured  along  central  Asian  paths  from  LTS  to  KNET. 

We  define  total  Pn  attenuation  as  the  combined  effect  of  geometrical  spreading  and 
finite  Q  and  calculate  it  at  a  reference  distance  of  1200  km  for  various  continental 
regions.  In  the  stable  regions  of  central  Asia,  Scandinavia  and  the  Canadian  shield  where 
Pn  velocity  is  higher  than  8.0  km/s,  Pn  attenuation  is  similarly  low.  Among  regions  in 
and  around  the  TP,  Pn  attenuation  is  the  strongest  under  north  central  Tibet,  and  decreases 
toward  the  south  and  east.  In  the  easternmost  TP,  Pn  attenuation  becomes  as  low  as  that  in 
the  stable  regions.  The  spatial  variation  of  Pn  attenuation  correlates  inversely  with  that  of 
Pn  velocity.  A  probable  cause  of  this  phenomenon  is  that  Pn  attenuation  is  directly 
affected  by  the  intrinsic  Qp  in  the  mantle  lid,  where  a  variation  in  temperature  causes 
both  Qp  and  vp  to  vary.  It  is  also  possible  that  regional  variations  of  velocity  structure  in 
the  vicinity  of  the  Moho,  such  as  variations  of  lid  gradient  and/or  density  and  aspect 
ratios  of  the  3D  scatterers,  also  play  a  role  in  causing  regionally  variable  Pn  attenuation. 
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Appendix 

Potential  Radiation  Pattern  from  the  1999  Earthquake 

In  theory,  Pn  spectra  from  earthquake  99.030  may  be  affected  by  a  non-isotropic 
source  radiation  pattern  (R(6 >)  in  equation  (2)).  In  this  appendix  we  demonstrate  that  the 
effect  of  this  pattern  is  small  and,  when  ignored,  should  not  have  caused  significant  errors 
in  Pn  Q  estimates.  First  we  note  that  this  study  uses  finite  Pn  wavetrains,  which  are  less 
vulnerable  to  the  radiation  pattern  than  are  the  first  Pn  pulses  (Zhao  and  Ebel,  1993).  To 
further  explore  the  level  of  the  radiation  pattern  from  event  99.030,  we  calculate  synthetic 
Pn  wavetrains  using  (1)  1-D  velocity  model  "Ml"  by  Roecker  et  al.  (1993),  and  (2)  the 
centroid  moment-tensor  (CMT)  solutions  by  Harvard  and  U.S.  Geological  Survey;  both 
solutions  characterizing  the  source  as  an  obliquely  reverse  rupture.  These  synthetics  are 
then  used  to  predict  Pn  radiation  pattern  to  the  directions  of  KNET  and  INDEPTH  III  sta¬ 
tions  (roughly  west  and  south,  respectively).  The  radiation  predicted  using  the  Harvard 
solution  is  about  30%  higher  to  KNET  than  to  INDEPTH  III  stations.  An  almost  exactly 
opposite  prediction  is  obtained  using  the  USGS  solution.  Perturbing  the  ID  velocity 
model  and  focal  solutions  cause  variations  of  the  predictions,  but  the  range  of  the  varia¬ 
tions  is  small,  presumably  because  neither  station  group  is  near  the  Pn  nodal  planes. 
Variations  in  radiation  pattern  of  such  level  (±30%  when  comparing  the  KNET  paths  to 
the  INDEPTH  paths)  are  less  than  the  observed  spectral  variations  within  each  station 
group  (Figure  4).  We  can  estimate  the  relative  error  in  Pn  Q{)  measurement  caused  by  an 
error  Sx  arising  from  an  unaccounted  R(9)  of  1.3,  using  equation  (A8)  of  Xie  et  al. 
(2004).  Assuming  <7x=0.3,  V=7.6  km/s,  /= 1  Hz,  and  A=1200  km,  we  obtain  SQ0/Q0  of 
16%  and  11%  when  the  true  Q0  is  260  and  175,  respectively.  These  relative  errors  in  Q0 
are  rather  small  and  acceptable.  In  the  true  Earth  structure,  3D  scattering  should  occur  to 
further  smooth-out  any  radiation  pattern  at  the  high-frequencies. 

Next  we  empirically  examine  whether  by  ignoring  the  effect  of  non-isotropic  R(6 >), 
our  estimates  of  Q 0  and  r]  to  the  INDEPTH  III  stations  are  significantly  biased  (see  sec¬ 
tion  "Estimates  of  Apparent  Pn  Q  to  Indepth  Stations").  Because  group  III-l  of  stations 
are  located  in  southern  Tibet  and  are  close  to  the  INDEPTH  II  stations  (here  after  referred 
to  as  "group  II";  see  Figures  1  and  5),  the  Pn  Q0  and  ij  from  the  October  7,  1994  explo¬ 
sion  (which  should  not  have  a  non-isotropic  R(9))  to  station  group  II  should  be  similar  to 
those  from  earthquake  99.030  to  station  group  III-1.  If  there  is  a  significantly  non- 
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isotropic  R(d)  by  event  99.030,  then  ignoring  it  should  result  in  a  biased  set  of  Q(]  and  rj 
estimates  over  paths  to  group  III- 1 .  The  bias  should  have  masked  the  similarity  between 
<20  and  7  values  to  group  III- 1  and  those  to  group  II.  In  the  main  text  we  ignored  R(6 >) 
and  estimated  average  Q0  and  r]  values  to  station  group  III- 1  to  be  255±48  and  0.1±0.1, 
respectively.  These  values  are  close  to  those  of  253±53  and  0.0±0.1,  estimated  for  paths 
to  station  group  II.  This  preserved  similarity  empirically  demonstrates  that  a  potential 
bias  in  the  the  estimated  Q0  and  7  values,  caused  by  a  non-isotropic  R{6)  from  event 
99.030,  is  small.  As  a  further  precaution,  we  repeated  the  procedure  of  estimating  Q0  and 
r 1  values  using  the  southern  half  (8)  stations  of  the  15-station  group  III- 1 .  These  8  stations 
are  even  closer  to  the  INDEPTH  II  stations  as  they  tend  to  overlap  (Figures  1  and  5).  We 
obtained  average  Q0  and  77  estimates  of  258±42  and  0.0±0.1,  respectively.  These  values 
are  similar  to  those  obtained  using  all  15  III- 1  stations  and  show  that  our  empirical  infer¬ 
ence  on  small  R(6)  is  robust. 
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Figure  Captions 


Figure  1.  Map  of  study  area.  See  inset  annotation  for  symbols  representing  tectonic 
boundaries  and  seismic  stations.  The  Sino-U.S.  (91-92)  and  INDEPTH  II  and  III  stations 
were  temporarily  deployed  during  three  PASSCAL  experiments  in  Tibet.  The  Khyrghis- 
tan  Network  (KNET)  is  a  permanent  network.  Regional  earthquakes  used  in  this  study  are 
shown  as  open  circles.  Earthquakes  whose  ID  numbers  start  with  91  and  92  occurred  in 
years  1991  and  1992  (see  Table  2  of  McNamara  et  al.,  1995  for  the  event  parameters). 
Earthquake  99.030  is  an  Mw= 5.5  (mb= 5.9)  event  that  occurred  on  January  30,  1999,  at 
3:51:5.4  UT.  Stars  are  the  June  10  and  October  7,  1994  Lop  Nor  explosions  (Table  1  of 
Xie  and  Patton,  1999).  White,  medium  and  dark  grays  represent  elevations  of  greater 
than  4000  M,  between  2000-4000  M,  and  less  than  2000  M,  respectively.  The  abbrevia¬ 
tions  are:  MTB  -  Main  Thrust  Belt,  IYS  -  Indus  Yalong  Suture;  BNS  -  Banggong-Nujiang 
Suture;  JRS  -  Jinsha  River  Suture,  KF  -  Kunlun  Fault,  ATF  -  Altyn  Tagh  Fault,  LTS  -  Lop 
Nor  Test  Site. 

Figure  2.  Vertical  component  record  section  from  the  INDEPTH  II  stations  recording 
the  October  7,  1994  explosion  (mb= 5.9).  Predicted  arrival  times  corresponding  to  approx¬ 
imate  Pn  and  Lg  group  velocities  of  7.6  and  3.5  km/s  are  marked  by  gray  lines.  Lg  is 
blocked.  See  Figure  3  of  Xie  and  Patton  (1999)  for  the  record  section  from  the  KNET. 

Figure  3.  Time  series  containing  the  Pn  wavetrain  and  the  preceding  noise  (top),  and 
Fourier  spectra  (bottom)  from  INDEPTH  II  station  BB15  recording  the  October  7  explo¬ 
sion.  For  clarity,  only  two  of  the  five  signal  windows  (the  first  and  third  windows)  are 
shown  by  the  dashed  curves  on  the  top.  Upper  and  lower  black  curves  in  the  bottom  are 
spectra  from  individual  signal  and  noise  windows.  The  spectrum  in  gray  is  the  average  of 
the  Pn  signal  spectra. 

Figure  4.  Pn  spectra  from  explosions  and  earthquake  in  the  Lop  Nor  region  of  eastern 

1  o 

Tarim  Basin.  Spectra  are  scaled  by  A  '  to  reduce  the  effect  of  geometrical  spreading. 
Each  row  is  from  one  event  as  indicated  at  lower  left,  with  "EX"  and  "EQ"  indicating 
explosions  and  earthquake  (see  Figure  1).  The  INDEPTH  III  stations  are  divided  into 
two  groups  (III- 1  in  row  3  and  III-2  in  row  4),  located  to  the  south  and  north  of  the  BNS 
(Figures  1  and  5).  Left  columns  show  individual  spectra  from  KNET  and  INDEPTH  net¬ 
works  as  black  and  gray  curves,  with  the  network  name  and  number  of  recording  stations 
written  nearby.  Right  columns  show  the  network  averaged  spectra  and  the  associated 
sample  standard  deviations.  The  average  distances  are  marked. 

Figure  5.  Map  showing  paths  used  in  this  study  and  the  estimated  Pn  Q(}  and  //  values. 
The  events  and  stations  are  the  same  as  those  plotted  in  Figure  1.  Black  paths  are  from 
earthquake  99.030  to  station  group  III-2,  sampling  the  northern  Tibetan  Plateau.  Light 
gray  paths  are  those  from  the  1994  explosions  to  INDEPTH  II  stations,  or  from  earth¬ 
quake  99.030  to  station  group  III- 1 ,  sampling  a  mixture  of  northern  and  southern  Tibetan 
Plateau.  White  paths  are  from  explosions  and  earthquake  99.030  to  KNET.  Dashed  lines 
plotted  in  darker  gray  are  paths  from  the  regional  earthquakes  to  the  1991-1992  PASS¬ 
CAL  stations,  used  in  the  two-station  Pn  Q  measurements.  Three  profiles  are  formed: 
SANG-TUNL,  WNDO-TUNL  which  is  a  subprofile  of  SANG-TUNL,  and 
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SANG/GANZ-MAQI  which  is  merged  from  two  sub-profiles  SANG-MAQI  and  GANZ- 
MAQI.  The  end  stations  along  these  profiles  (SANG,  TUNL,  WNDO,  SANG,  GANZ 
and  MAQI)  are  written.  See  McNamara  et  al.  (1995)  for  other  station  names.  The  Pn  Q0 
and  q  values  are  written  near  the  paths  over  which  they  are  estimated  (values  for  paths  in 
gray  are  the  average  to  station  groups  II  and  III- 1);  oblique  numbers  are  those  obtained 
using  the  two-station  method. 

Figure  6.  Network  averaged  Pn  spectra  (fluctuating  curves),  and  the  model  spectra 
(smooth  curves)  corresponding  to  rows  1,  3  and  4  of  Figure  4.  Model  spectra  are  con¬ 
structed  using  the  MMM  or  Brune  source  model,  the  estimated  source  M0,  fc  and  path- 
averaged  <2o  and  q  values.  Similarly  to  Figure  4,  black  and  gray  are  used  for  spectra  from 
KNET  and  INDEPTH  networks.  The  network  name  and  path-averaged  Q0  and  7  values 
are  written  near  the  curves.  Note  that  only  8  KNET  stations  recorded  the  October  7 
explosion,  so  the  average  Q0  is  381,  rather  than  364  obtained  using  all  10  KNET  stations 
(Xie  and  Patton,  1999). 

Figure  7.  Stacked  spectral  ratios  (SSRs)  calculated  along  three  profiles  in  eastern  Tibet 
(dots)  and  the  best-fitting  models  of  Q0  and  q  (straight  lines).  Profile  names  and  numbers 
of  two-station  pairs  used  are  written  on  the  top  of  each  panel  (see  Figure  5  for  profile 
locations).  The  frequency  bands  used,  and  estimated  Q0  and  q  values  are  written  inside 
the  panels. 

Figure  8.  Pn  attenuation  models  calculated  for  various  regions  from  published  observa¬ 
tions,  at  a  reference  distance  of  1200  km.  Each  model  is  plotted  in  the  frequency  range  of 
observation.  Model  "Canadian  S."  is  for  the  Canadian  shield,  calculated  using  the  Q0  and 
q  values  and  the  G.S.T.  by  Zhu  et  al.  (1991).  Model  "Scandinavia"  is  from  Sereno  et  al. 
(1988).  Model  "C.  Asia"  is  for  central  Asia  (Xie  and  Patton,  1999;  this  study).  Models 
"C.  TP."  and  "N.  TP."  are  obtained  from  this  study  for  paths  running  from  LTS  to  stations 
in  south-central  Tibet  and  north-central  Tibet  (gray  and  black  paths  in  Figure  5),  respec¬ 
tively.  Values  of  Q0  and  q  of  the  models  obtained  in  this  study  are  written.  Models  for 
east  Tibet  (along  profile  SANG-TUNL  in  Figure  5)  and  easternmost  Tibet  (along  profile 
SANG/GANZ-MAQI),  both  obtained  in  this  study  are  not  plotted  because  they  are  simi¬ 
lar  to  Models  C.  TP  and  C.  Asia  (see  section  "Apparent  Pn  Q  along  Passcal  Profiles  in 
Eastern  Tibet"),  respectively. 

Figure  9.  Stacked  spectral  ratios  (SSRs)  calculated  along  profile  SANG-TUNL  using  a 
G.S.T.  correction  of  m=l.  07  +  0. 035/  (in  gray),  versus  SSRs  with  a  G.S.T.  correction  of 
m=1.3  (in  black;  the  same  as  in  the  top  panel  of  Figure  7).  The  best  fit  Q{)  and  q  models 
and  resulting  residuals  are  indicated.  The  residuals  differ  from  each  other  by  less  than 
0.2%,  indicating  a  virtually  complete  trade-off  between  G.S.T.  and  Q0  and  q  for  the  data 
used  in  this  study.  The  former  G.S.T.  is  estimated  for  the  Canadian  shield  by  Zhu  et  al 
(1991)  and  is  used  here  to  illustrate  the  parameter  trade-off. 
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Xie  2006,  Fig.  1. 


Xie  2006,  Fig.  5. 


Xie  2006,  Fig.  6. 
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Xie  2006,  Fig.  7. 
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Eurasia  from  Rayleigh-wave  Phase  and  Amplitude  Data 

James  B.  Gaherty  \  Po  Chen  \  Li  Zhao2,  and  Tae-Kyung  Hong1’3 

1  Lamont-Doherty  Earth  Observatory  of  Columbia  University 
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Abstract 

We  have  mapped  lateral  variations  of  crustal  and  upper-mantle  shear- wave 
velocity  ft  and  shear  attenuation  (1  IQf)  in  eastern  Eurasia  using  a  new  joint  inversion 
procedure  based  on  three-dimensional  partial-derivative  (Frechet)  kernels.  In  our 
analysis,  we  utilized  Rayleigh  waves  with  periods  longer  than  —15  s,  derived  from 
regional  and  near-teleseismic  earthquakes.  We  collected  and  analyzed  859  vertical- 
component  seismograms,  derived  from  over  133  events  with  Mw  >5.5  that  occurred 
within  the  eastern  Eurasia  region  (10°-60°N,  70°-140°E)  between  2000-2006.  The  data 
were  recorded  on  71  broadband  seismic  stations  within  eastern  Eurasia  (10°-60°N,  80°- 
140°E).  The  seismograms  were  evaluated  for  signal-to-noise  within  the  band  of  interest 
(10  mHz  -  80  mHz)  by  comparing  the  records  to  complete  normal-mode  synthetic 
seismograms  for  radially  symmetric  and  transversely  anisotropic  reference  models.  We 
then  made  15,462  mutually  consistent  measurements  of  frequency-dependent  traveltime 
and  amplitude  anomalies  for  surface  waves  by  cross-correlation  between  records  and 
synthetics.  We  used  coupled  normal-mode  summation  to  compute  the  three-dimensional 
Frechet  kernels  of  these  traveltime  and  amplitude  anomalies  to  both  the  elastic  shear- 

wave  speeds  fl  and  shear- wave  quality  factor  0W  values.  Our  starting  model  for  Qu  was 
extracted  from  a  recent  global  model  of  upper-mantle  attenuation.  The  frequency- 
dependent  traveltime  and  amplitude  anomalies  were  jointly  inverted  for  three- 
dimensional  models  of  both  elastic  and  anelastic  structures.  This  coupled  inversion 
allowed  us  to  account  for  the  effects  on  amplitudes  from  elastic  heterogeneities  such  as 
scattering,  focusing  and  defocusing,  as  well  as  the  effects  of  anelastic  dispersion  on  the 
traveltimes.  The  regional  models  we  have  obtained  provide  lateral  resolution  of  about 
400  km  in  the  crust  and  upper  mantle.  In  the  crust,  high  values  of  Qu  are  concentrated  in 
eastern  China,  extending  southward  from  the  north  China  craton,  while  a  clear  low-0 
anomaly  in  the  crust  beneath  northern  Tibet.  In  the  upper  mantle,  the  model  is 
characterized  by  a  lithosphere  with  very  high  0  values  ranging  from  250-400,  underlain 
by  a  low-0  asthenosphere  with  0^  of  60-100.  The  lateral  variation  of  the  upper-mantle  Qu 
perturbation  appears  to  correlate  with  the  lateral  variation  in  the  inverted  shear-velocity 
perturbation. 

Introduction 

Accurate  event  detection,  identification,  and  characterization  are  critical 
components  of  the  US  nuclear  monitoring  program.  A  necessity  for  improving  these 
estimates  is  the  more  accurate  estimation  of  seismic  attenuation  rate  (or  its  inverse,  0), 
which  controls  the  amplitudes  of  seismic  phases  upon  which  the  assessments  are  based. 
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Q  is  known  to  vary  widely  in  continental  regions,  making  reliable  estimates  of  it  critically 
important,  but  also  quite  difficult.  Q  is  also  frequency  dependent,  which  makes  it  difficult 
to  determine  whether  Q  estimates  derived  from  different  phases  with  different  frequency 
characteristics  are  a  result  of  intrinsic  frequency  dependence,  or  differences  in  structural 
sampling.  Finally,  Q  estimates  are  derived  from  seismic  amplitudes,  which  are  also 
strongly  dependent  on  elastic  effects  such  as  wavefield  focusing  and  scattering.  As  a 
result,  robust  tomographic  (3D)  quantification  of  Q  on  regional  scales  remains  in  its 
infancy. 

We  address  this  problem  by  developing  and  applying  a  comprehensive  waveform 
analysis  and  inversion  scheme  to  simultaneously  estimate  3D  variations  in  shear  velocity 
(/?)  and  shear-wave  attenuation  (or  more  precisely,  its  inverse,  the  shear-wave  quality 
factor  Qlt)  from  Rayleigh  waves  recorded  in  eastern  Eurasia.  The  analysis  technique 
characterizes  phase  and  amplitude  residuals  between  recorded  and  synthetic  Rayleigh 
waves,  and  inverts  them  using  new  three-dimensional  sensitivity  (Frechet)  kernels  for 
structural  perturbations.  The  analysis  spans  a  frequency  band  from  0.01-0.05  Hz,  with  the 
low  frequency  estimates  providing  a  stable  long-wavelength  framework  from  which  we 
bootstrap  to  higher  frequency.  We  focus  on  Rayleigh  waves  because  they  provide  self- 
consistent  spatial  and  depth  constraints  on  Qu  in  the  crust  and  upper  mantle  across  a 
broad  frequency  band.  In  Eurasia,  Rayleigh  waves  from  regional  and  near-teleseismic 
events  give  good  spatial  coverage,  and  can  quantify  large-scale  regional  trends  and  the 
effects  of  long-wavelength  scattering.  The  resulting  model  spans  all  of  China,  the  Korean 
peninsula,  and  much  of  Southeast  Asia,  and  can  be  used  to  place  localized  high- 
frequency  observations  (e.g.  Lg,  Pg,  S„,  P„)  into  a  larger  regional  context. 

Data  Collection 

Our  data  set  consists  of  broadband  Rayleigh  waves  that  traverse  Eastern  Eurasia. 
The  seismic  records  were  recorded  on  instruments  from  the  IRIS  GSN,  cooperating 
networks  such  as  the  CDSN,  and  several  IRIS  PASSCAL  experiments  from  Asia.  All 
data  were  collected  from  IRIS  DMC.  The  stations  used  are  distributed  in  the  region  from 
10°N  -  60°N  and  80°E  -  140°E.  We  used  stations  recording  either  an  LHZ  or  BHZ 
component,  with  a  total  of  71  stations.  Events  occurred  in  the  region  between  10°N  - 
60°N  and  70°E  -  140°E  in  the  time  period  between  2000  and  2006.  The  magnitudes  are 
greater  than  Mw  5.5.  The  initial  dataset  was  cleaned  to  remove  data  with  a  signal-to-noise 
level  of  less  than  approximately  2  in  the  central  portion  of  the  frequency  band  (0.02-0. 1 
Hz),  as  well  as  data  clearly  contaminated  by  glitches  or  other  problems.  The  final  number 
of  records  used  is  859  from  133  events.  The  final  event  and  station  distribution  as  well  as 
the  spatial  path  coverage  associated  with  this  data  is  shown  in  Figure  1. 

Methodology 

We  analyzed  the  selected  seismograms  using  the  generalized  seismological  data 
functionals  (GSDF)  method  (Gee  and  Jordan,  1992;  Gaherty  et  ah,  1996),  where  the 
phase-delay  times  3tfcon)  and  differential  amplitudes  8tf(on)  of  observed  Rayleigh  waves 
were  measured  relative  to  complete  normal-mode  synthetic  seismograms  as  a  function  of 
frequency  a>n  across  the  0.01-0.05  Hz  band.  This  methodology  utilizes  synthetic 


35 


seismograms  of  a  target  wavegroup  (in  this  case  the  fundamental-mode  Rayleigh  waves) 
to  estimate  phase  delays  and  relative  amplitudes  of  the  observed  arrival  relative  to  the 
synthetic  as  a  function  of  frequency.  An  example  of  the  GSDF  processing  is  shown  in 
Figure  2.  To  make  the  GSDF  measurements,  the  synthetic  wavegroup,  called  an 
“isolation  filter”,  is  cross-correlated  with  both  the  data  and  a  complete  synthetic 
seismogram.  The  resulting  cross-correlagrams  are  windowed  and  narrow-band 
(Gaussian)  filtered  at  discrete  frequency  intervals  between  0.01  and  0.05  Hz,  and  the  peak 
phase  and  amplitude  of  each  correlagram  is  estimated  at  each  frequency.  Referencing  the 
observations  to  a  synthetic/isolation-filter  cross-correlation  in  this  fashion  allows  us  to 
account  for  interference  from  unmodeled  wavegroups  and  minimizes  bias  associated  with 
windowing  and  filtering.  The  synthetics  are  complete  from  0-50  mHz  (all  periods  longer 
than  12.5  s),  and  are  calculated  for  a  reference  models  (Figure  3)  constructed  from  a 
composite  of  IASP91  for  mantle  shear  and  compressional  velocities;  Crust  2.0  for  local 
crustal  velocity  perturbations  (Laske  et  ah,  2001);  a  mantle  (9/,  profile  for  the  China 
region  extracted  from  a  recent  global  model  (Dalton  and  Ekstrom,  2006);  and  a  crustal  QLl 
profile  derived  from  analysis  of  regional  phases  (Jemberie  and  Mitchell,  2004). 

These  observations  characterize  the  observed  wavefield  characteristics  in  a 
manner  similar  to  traditional  estimates  of  phase  velocity  and  spectral  amplitudes.  They 
differ  from  the  traditional  estimates  in  that  they  are  relative  to  a  reference  model  (similar 
to  travel-time  residuals),  and  thus  can  be  directly  inverted  for  improved  velocity  and  Q 
structure  using  a  linearized  Gaussian-Bayesian  approach  (e.g.  Gaherty  et  ah,  1996).  In  the 
frequency  domain,  we  can  map  the  synthetic  waveform  u^oj)  into  the  observed 
wavefonn  u^co)  using  two  frequency-dependent,  time-like  quantities  <5YP( co)  and  Srq(co) 

Uj(co)  =  u;(co)  exp  {/'  co[  8  rp  ( co )  +  /dVq(^)]|.  (1) 

The  GSDF  measurements  of  8tx(oon)  (x  =  p,  q)  are  weighted  averages  of  8rx  over  a  narrow 
frequency-band  centered  on  con.  As  shown  in  Chen  et  al.  (2006),  these  frequency- 
dependent  GSDF  measurements  can  completely  capture  the  waveform  information  if  the 
sampling  in  frequency  domain  is  dense  enough.  And  the  GSDF  measurements  are  well 
suited  for  the  tomographic  inverse  problem.  In  particular,  their  linearization  is  based  on 
the  Rytov  approximation,  which  is  valid  for  large  accumulative  phase-shifts  as  long  as 
the  phase  perturbation  per  wavelength  is  small  (Chernov  1960;  Snieder  &  Lomax  1996). 
This  is  much  less  restrictive  than  the  Born  approximation,  which  requires  small 
accumulative  phase-shifts. 

The  sensitivity  (Frechet)  kernels  for  our  GSDF  measurements  with  respect  to 
shear- wave  velocity  /3  and  attenuation  QM  were  computed  using  the  “full  wave”  approach 
of  Zhao  &  Jordan  (2006).  First,  we  linearize  the  relation  between  GSDF  measurements 
and  waveform  perturbation 

SdZ=)dtJZ(t)Su-(xr,t),  (2) 
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where  ds-n  as  the  /7-th  GSDF  measurement  made  on  the  7-th  component  of  the 
seismogram,  which  is  generated  by  source  s  and  recorded  at  receiver  r.  Exact  expressions 
of  the  seismogram  functional  Jsd(t)  for  GSDF  measurements  are  give  in  Chen  (2005). 
Second,  we  linearize  the  relation  between  waveform  perturbation  and  the  perturbation  to 
the  density  8p(x)  and  elastic/anelastic  modulus  tensor  8cjklm  (x)  using  the  Born 
approximation 


Su.(xr,t)  =  -\  dV(x)  J  d rj (xr,t~  z;x) 8] w' (x, z) 8p(x) 

+ Tjdk GtJ  (xr,  t  -  z;x)  d,  um  (x,  z)8cJklm  (*)]  > 
jklm  J 


(3) 


The  exact  Frechet  kernels  for  GSDF  measurements  can  therefore  be  expressed  as 
KPdM  =  -  j  dt\  d xJl ( t)X Gn (\,t-z;\r )82r u) ( x,  r) ,  (4) 

j 

K'jddx)  =  -  j  dl  \  dzdli^Gfi/xJ  -  T-,x  r)G,um(x,z),  (5) 

which  generally  involves  the  convolution  between  gradients  of  the  forward  earthquake 
wavefields  um(\J)  with  the  gradient  of  the  transposed  “receiver  Green  tensor”  (RGT) 
GC(x,t;xr)  (Zhao  et  al.  2006).  Because  both  the  forward  earthquake  wavefields  and  the 

RGTs  were  computed  using  the  normal-mode  theory,  which  provides  the  exact  solution 
to  the  wave  equation  in  a  ID  reference  model,  the  Frechet  kernels  constructed  this  way 
can  account  for  the  full  physics  of  wave  propagation  including  multipathing  and  non- 
geometrical  arrivals  (Zhao  &  Jordan  2006). 

When  intrinsic  attenuation  is  taken  into  account,  the  (real)  elasticity  tensor  is 
replaced  by  a  complex  tensor: 

cMm = |v(®)[i + ;o»]  -  + iQ-;m\sA 

+  //(<»)[  1  +  iQj  (  tu)  ]  ( Sj[  8 km  +  Sjm8kl)  +  yjklm, 

where  QK  and  Qu  are  the  quality  factors  for  the  incompressibility  k  and  the  shear 
modulus  /j,  respectively.  In  (6),  the  elasticity  tensor  has  been  decomposed  into  a  purely 
isotropic  part  expressed  in  tenns  of  k  and  //  and  a  purely  anisotropic  part  yjkIm ,  and  the 

intrinsic  attenuation  of  only  the  purely  isotropic  part  is  considered  (Dahlen  &  Tromp 
1998).  Substituting  the  complex  tensor  in  equation  (6)  into  equation  (3),  we  can  obtain  a 
linearized  relation  between  the  waveform  perturbation  Susi(xr,t)  and  the  perturbations 
8Qk  and  8Qn .  From  this  linearized  relation  and  equation  (2)  we  can  derive  the  expressions 
for  the  Frechet  kernels  of  the  GSDF  measurements  with  respect  to  the  shear-wave 
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velocity  [3  and  attenuation  factor  Q  .  Figure  4  shows  examples  of  Frechet  kernels  of 

phase-delay  time  8tv  with  respect  to  shear-wave  velocity  and  amplitude-reduction  time 
Stq  with  respect  to  shear-wave  velocity  and  Q  for  a  Rayleigh  wave  at  35  mHz. 


Results 

We  made  7,731  frequency-dependent  phase-delay  time  ( Stp )  and  7,731  amplitude 
anomaly  measurements  (<5/q)  for  the  Rayleigh  waves  on  the  selected  859  high-quality 
vertical-component  seismograms  at  9  frequencies  ranging  from  10  mHz  to  50  mHz.  The 
15,462  GSDF  measurements  were  then  inverted  jointly  for  both  shear-wave  velocity 
perturbation  S/3//3  and  attenuation  perturbation  (3QJ  using  the  full-wave  Frechet 
kernels  (Figure  4).  This  joint  inversion  for  both  shear  velocity  and  attenuation  allows  us 
to  account  for  the  effects  of  scattering  and  focusing/defocusing  effects  on  waveform 
amplitudes  cause  by  elastic  heterogeneities  using  the  amplitude  kernels  with  respect  to 
shear  velocity.  The  phase-delay  measurements  included  in  the  inversion  provide 
additional  constraints  on  shear  velocity.  In  this  model,  we  do  not  incorporate  the 
frequency-dependence  of  Qm;  the  derived  models  are  appropriate  for  a  center  frequency 
of  ~25  mHz,  and  assume  constant  Q  over  a  frequency  band  of  0.01-0.06  Hz. 

Our  modeling  volume  extends  from  70°  E  to  141°  E  in  longitude,  from  10°  N  to 
58°  N  in  latitude  and  from  0  km  to  280  km  in  depth.  We  discretized  the  modeling  volume 
into  2°  by  2°  blocks  with  a  vertical  dimension  varying  from  10  km  at  shallower  depth  to 
40  km  at  larger  depth.  The  model  parameters  are  assumed  to  be  constant  within  each 
block.  This  parameterization  results  in  21,600  model  parameters  considering  both  shear 
velocity  perturbation  5/3/ (3  and  attenuation  perturbation  dQ^J  0L1.  The  sensitivity  kernels, 
which  were  originally  computed  on  a  1 0  by  1 0  grid  with  vertical  grid  spacing  of  1  km, 
were  then  averaged  over  each  2°  block.  We  introduced  266  additional  model  parameters 
to  account  for  the  unknown  errors  in  source  origin  time  and  scalar  moment  for  the  133 
earthquakes  used  in  this  study.  The  total  number  of  model  parameters  used  in  our 
inversion  is  21,866. 

The  resulting  linear  system  was  then  scaled  using  a  data- weighting  matrix  Cd12 

and  concatenated  with  the  “roughing  operator”  Cm12 ,  which  is  a  linear  combination  of 
the  Laplacian  operator  (Constable  et  al.  1987;  Sambridge  1990;  Tarantola  2005)  and  the 
identity  operator,  i.e.  Cm12  =  0(\  +  TV2).  Using  this  definition  of  the  “roughing  operator”, 
we  impose  our  prior  infonnation  on  the  inversion;  i.e.,  we  assume  that  in  the  absence  of 
any  other  information,  the  model  perturbation  is  smooth  and  small.  We  note  that  the 
corresponding  “smoothing  operator”  Cjj2  can  be  generated  from  an  exponential 
correlation  function,  whose  correlation  length  is  proportional  to  A.  In  practice,  the 
Laplacian  operator  is  approximated  numerically  by  finite  differencing.  The  final  linear 
system  that  we  solve  is 


C,1/2A 

c5m  = 

T5 

<N 

l) 

d 

a 

^1— 1/2  I 
-  '“'m 

0 
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where  matrix  A  contains  the  averaged  sensitivity  kernels  for  each  GSDF  measurement, 
8m  is  the  model  parameter  vector  and  d  is  the  data  vector.  The  linear  system  specified  by 
equation  (7)  was  solved  using  a  parallelized  LSQR  subroutine  from  the  PETSc  parallel 
scientific  computing  library  (Balay  et  al.  2004).  The  LSQR  method  (Paige  &  Saunders 
1982),  which  is  a  type  of  conjugate-gradient  method,  provides  a  very  efficient  means  for 
solving  large-scale  sparse  linear  systems. 

We  have  carried  out  a  series  of  inversions  to  construct  the  trade-off  curves 
between  variance  reduction  and  model  smoothness  (Figure  5).  Our  preferred  model 
EAF3D  locates  around  the  corners  of  the  trade-off  curves.  Figure  6  shows  this  preferred 
model  at  12  different  depths.  This  model  provides  a  variance  reduction  of  about  52%  and 
is  relatively  smooth. 

An  important  indicator  of  the  resolution  of  our  inversion  is  the  spatial  coverage  of 
all  the  Frechet  kernels  used  in  the  inversion.  In  Figure  7,  we  have  plotted  the  summation 
of  the  absolute  values  of  all  the  kernels  within  each  block  of  our  modeling  volume.  The 
best  kernel  coverage  is  in  the  crust  and  upper  mantle;  the  coverage  is  still  decent  down  to 
a  depth  of  about  170  km. 

Discussion 

In  general,  the  inverted  model  perturbation  EAF3D  is  geologically  reasonable. 

For  shear-velocity  perturbation,  the  amplitude  ranges  from  -5%  to  5%.  Spatially,  the 
velocity  model  shows  a  broad  region  of  low  shear- wave  speeds  in  central  and  western 
China  that  extend  through  the  crust  and  upper  mantle,  with  higher  velocities  to  the  east. 
These  velocity  anomalies  are  derived  primarily  to  account  for  their  effect  on  wave 
amplitude  (and  thus  Q),  and  thus  are  not  the  focus  of  this  study.  For  Qu  perturbation,  the 
amplitude  ranges  from  -40%  to  40%.  In  the  crust,  high  values  of  Qu  are  concentrated  in 
eastern  China,  extending  southward  from  the  north  China  craton.  This  result  is  in  good 
agreement  with  variations  observed  by  Jemborie  and  Mitchell  (2004)  from  Rayleigh 
waves,  and  also  generally  agrees  with  the  result  of  Xie  et  al.,  (2005)  derived  from  high 
frequency  Lg  observations.  A  clear  low-0  anomaly  in  the  crust  beneath  northern  Tibet 
is  also  in  good  agreement  with  previous  studies  (Jemborie  and  Mitchell,  2004;  Xie  et  al., 
2005),  and  the  anomaly  persists  well  into  the  upper  mantle  (Fig  6b). 

A  general  feature  of  our  inverted  QM  perturbation  is  that  at  shallower  depth  in  the 
lithosphere,  the  average  perturbation  is  positive;  while  at  larger  depth  in  the 
asthenosphere,  the  average  perturbation  is  negative.  This  inverted  QM perturbation  seems 
to  enhance  the  depth  gradient  that  already  exists  in  the  ID  reference  model  (Figure  3).  If 
we  apply  the  inverted  QM  perturbation  to  the  ID  reference  model,  we  obtain  an  updated 
3D  model  that  characterizes  both  the  lateral  and  depth  variation  of  the  absolute  value  of 
(^throughout  the  region  (Figure  8).  In  the  crust,  varies  from  a  minimum  of  75  to  a 
maximum  near  200.  This  is  similar  to  the  range  observed  by  Jemborie  and  Mitchell 
(2004),  although  in  detail  their  minimum  value  beneath  Tibet  is  slightly  slower  ( — 4 0 ) , 
and  our  maximum  values  beneath  eastern  China  are  slightly  higher  than  theirs.  These 
differences  are  likely  due  to  resolution  differences  in  the  models.  Jemborie  and  Mitchell 
likely  contain  better  resolution  in  the  upper  crust,  while  our  explicit  inclusion  of  Q 
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variations  in  the  mantle  may  result  in  improved  Q M  estimates  in  the  lower  crust.  Further 
tests  of  these  models  are  underway  to  evaluate  this  issue. 

In  the  upper  mantle,  the  model  is  characterized  by  very  high  Q  values  ranging 
from  250-400.  This  is  largely  driven  by  the  high -Q  region  in  eastern  China.  The  result  is 
robust,  in  that  the  starting  model  had  a  relatively  high -Q  lid  (Figure  3),  but  the  relative 
amplitude  residuals  for  the  Rayleigh-wave  frequencies  sensitive  to  lid  structure  clearly 
required  even  higher  Q  on  average.  As  a  result  of  model  smoothness  constraints,  the 
absolute  value  of  Q  in  the  low  Q  region  of  western  China  is  likely  overestimated.  Below 
the  lid,  the  Q  drops  to  very  low  values  in  the  range  of  60-100.  Again,  on  average  these 
values  seem  to  be  quite  robust,  as  the  amplitudes  of  the  long-period  Rayleigh  waves  were 
persistently  low,  despite  the  low-Q  asthenosphere  in  the  starting  model.  This  sharp 
transition  from  a  high-Q  lid  to  a  low-Q  asthenosphere  appears  to  persist  throughout  most 
of  China  at  depths  of  135  km  and  greater. 

The  lateral  variation  of  the  upper-mantle  ^perturbation  seems  to  correlate  with 
the  lateral  variation  in  the  inverted  shear-velocity  perturbation.  This  result  is  consistent 
with  the  conclusion  reached  by  Xie  (2006)  based  on  spectral  analysis  on  Pn  arrivals  from 
1  Hz  to  10  Hz.  In  Xie  (2006),  this  lateral  variation  of  attenuation  is  interpreted  as  a  result 
of  temperature  variation  in  the  mantle  lid  under  Tibet,  i.e.  the  temperature  in  the  mantle 
lid  under  northern  Tibet  is  abnormally  high  and  causes  partial  melting  while  it  is 
abnormally  low  under  southern  and  eastern  Tibet  due  to  the  underthrusting  Indian 
lithosphere. 

We  emphasize  that  these  models  are  preliminary.  We  are  currently  updating  the 
models  by  incorporating  two  improvements: 

1 .  We  are  incorporating  the  effects  of  variable  topography  of  major  discontinuities 
such  as  the  Moho.  Explicit  expressions  for  the  Frechet  kernels  of  the  waveform 
with  respect  to  discontinuity  topography  are  given  in  Dahlen  (2005).  We  can 
incorporate  this  information  into  our  tomographic  inversion  either  as  explicit 
constraints  in  our  model  parameterization  or  correct  for  the  effects  of  variable 
topography  on  our  GSDF  measurements  using  the  exact  Frechet  kernels.  We  are 
testing  both  approaches 

2.  We  will  increase  the  upper  bound  of  our  inverted  frequency  range  from  50  mHz 
to  80  mHz.  Currently,  our  synthetics  and  measurements  extend  to  80  mHz,  but  the 
high-frequency  observations  are  still  being  evaluated  for  stability.  They  have  not 
been  included  in  the  models  presented  here.  Once  these  data  have  been  evaluated 
(specifically,  corrected  for  the  effects  of  cycle  skipping),  they  will  be  included  in 
the  inversion. 

Conclusions 

We  have  developed  a  comprehensive  waveform  analysis  and  inversion  scheme  to 
simultaneously  estimate  3D  variations  in  shear  velocity  (//)  and  shear-wave  quality  factor 
Q/j,  and  applied  this  scheme  to  a  large  Rayleigh  waves  dataset  sampling  eastern  Eurasia. 
The  analysis  technique  characterizes  phase  and  amplitude  residuals  between  recorded  and 
synthetic  Rayleigh  waves,  and  inverts  them  using  new  three-dimensional  sensitivity 


40 


(Frechet)  kernels  for  structural  perturbations.  The  joint  inversion  allows  us  to 
incorporate  the  effects  of  wavefield  focusing  and  scattering  on  wave  amplitudes, 
improving  our  ability  to  resolve  variations  in  attenuation.  The  analysis  spans  a  frequency 
band  from  0.01-0.05  Hz,  with  the  low  frequency  estimates  providing  a  stable  long- 
wavelength  framework  from  which  we  bootstrap  to  higher  frequency.  The  resulting 
model  provides  self-consistent  spatial  and  depth  constraints  on  Qu  in  the  crust  and  upper 
mantle,  with  spatial  resolution  of  several  hundred  kilometers.  Spatial  variations  of  (Qyy 
observed  in  the  crust  generally  agree  with  those  observed  from  higher-frequency  crustal 
phases.  Spatial  variations  in  the  mantle  Q M  correlate  with  shear  velocity,  with  high  0Lt 
lithosphere  observed  in  eastern  China,  and  a  low  Qyi  region  found  beneath  Tibet.  The 
average  attenuation  across  the  region  appears  to  be  highly  stratified,  with  a  very  high  Qyy 
lithosphere  underlain  by  very  low  Qyt  in  the  asthenosphere. 
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Figure  1.  Sour-station  path  coverage  for  our  current  inversion.  Harvard  CMT  solutions 
are  plotted  at  the  epicenters  of  the  earthquakes.  Stations  are  shown  as  blue  triangles.  Red 
lines  shows  source-station  paths  used  in  our  inversion. 
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Figure  2.  An  example  of  GSDF  data  processing,  (a)  observed  seismogram;  (b)  complete 
synthetic  seismogram,  the  isolation  filter  is  shown  in  red;  (c)  cross-correlagram  between 
the  isolation  filter  and  the  observed  seismogram;  (d)  cross-correlagram  between  the 
isolation  filter  and  the  complete  synthetic  seismogram;  (e)  frequency-dependent  phase- 
delay  measurements;  (f)  frequency-dependent  amplitude  anomaly  measurements. 
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Figure  4(a).  Partial  derivative  kernels  for  phase-delay  with  respect  to  shear  velocity,  for 
a  typical  Rayleigh  wave  observation,  (top  panel)  Broadband  (5-60  mHz)  Rayleigh 
waveform  selected  for  analysis  (dashed  lines  represent  window),  (middle  panel)  Cross- 
section  of  kernel  within  the  vertical  source-receiver  plane,  (bottom  panel)  Plan  view  of 
kernel  at  depth  20  km. 
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Figure  4(b).  Partial  derivative  kernels  for  relative  amplitude  with  respect  to  shear 
velocity,  for  the  example  Rayleigh-wave  observation  in  Figure  1(a).  (top  panel)  Cross- 
section  of  kernels  within  the  vertical  source-receiver  plane,  (bottom  panel)  Plan  view  of 
kernel  at  depth  20  km. 
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Figure  4(c).  Partial  derivative  kernels  for  relative  amplitude  with  respect  to  Q M,  for  the 
same  Rayleigh  wave  example  from  Figure  1(a).  (left)  Cross-sections  of  kernels  within 
the  vertical  source-receiver  plane,  as  well  as  perpendicular  to  the  path,  midway  between 
source  and  receiver,  (right)  Plan  view  of  kernels  at  depths  of  5,  20,  and  100  km.  In  all 
kernels,  color  scale  is  such  that  blue  indicates  regions  where  an  increase  in  will 
produce  an  increase  in  predicted  Rayleigh-wave  amplitude 
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Figure  5.  The  trade-off  relations  between  model  size  or  smoothness  and  variance 
reduction.  Variance  reduction  as  a  function  of  the  model  norm  damping  parameter  6  for 
shear  velocity  (a);  of  the  model  smoothness  damping  parameter  A  for  shear  velocity  (b); 
of  the  model  norm  damping  parameter  0  for  QM  (c);  of  the  model  smoothness  damping 
parameter  A  for  Qu. 
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Figure  6(a)  Shear-velocity  perturbation  (left)  and  attenuation  perturbation  (right)  at  5  km 
to  3 1  km  depths. 
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Figure  6(b)  Shear-velocity  perturbation  (left)  and  attenuation  perturbation  (right)  at  49 
km  to  136  km  depths. 
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Figure  6(c)  Shear-velocity  perturbation  (left)  and  attenuation  perturbation  (right)  at  206 
km  to  28 1  km  depths. 
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Figure  7  Kernel  coverage  at  12  different  depths.  First  and  third  columns:  coverage  of 
shear-velocity  kernels;  second  and  fourth  columns:  coverage  of  attenuation  kernels. 
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Figure  8(a)  The  updated  3D  shear-wave  attenuation  model  for  eastern  Eurasia  obtained 
by  applying  the  3D  perturbation  (Figure  6)  to  the  ID  reference  model,  for  depths  of  5-135 
km. 
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Figure  8(b)  The  updated  3D  shear-wave  attenuation  model  for  eastern  Eurasia  obtained 
by  applying  the  3D  perturbation  (Figure  6)  to  the  ID  reference  model,  for  depths  of  169- 
281  km. 
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